Signal processing

ABSTRACT

A particular computer-implemented method includes receiving sensed data from sensors of a sensor array, where data from each sensor is descriptive of waveform phenomena detected at the sensor. The method also includes determining an estimated spatial spectrum of the waveform phenomena based at least partially on the sensed data. The method further includes determining an estimated covariance matrix of the waveform phenomena based on the estimated spatial spectrum. The method includes determining adaptive beamforming weights using the estimated covariance matrix.

CLAIM OF PRIORITY

The present application claims priority to U.S. Provisional ApplicationNo. 61/408,456 filed Oct. 29, 2010, which is incorporated by referenceherein in its entirety.

FIELD OF THE DISCLOSURE

The present disclosure is generally related to signal processing.

BACKGROUND

Deployment of long aperture arrays using large numbers of sensorelements has been enabled by the evolution of sensor technology,telemetry technology, and digital processing technology. Such arrayshold promise for benefits in terms of theoretical processing gain andspatial resolution. Long aperture arrays have been employed across arange of applications, including acoustic and radar applications.Adaptive beamforming techniques have been applied to long aperturearrays to increase a system's ability to perform well in complicatedsignal and interference environments. For many adaptive beamformingtechniques, increasing the number of sensor elements of the array worksagainst a desire to adapt quickly. Additionally, multipath environmentsintroduce challenges that may cause substantial performance loss forcertain adaptive beamforming techniques.

Many approaches for covariance estimation for adaptive beamforming havebeen proposed. Many adaptive beamforming techniques use sensed data toestimate a covariance matrix and use the covariance matrix for spectralestimation. However, snapshot deficient operation and correlated signaland interference environments continue to be difficult to handle forcertain adaptive beamforming techniques.

SUMMARY

An adaptive beamforming technique is disclosed in which a spatialspectrum is estimated based on sensed data. The estimated spatialspectrum is used to estimate a covariance matrix. The estimatedcovariance matrix may be used for adaptive beamforming.

In a particular embodiment, a computer-implemented method includesreceiving sensed data from sensors of a sensor array, where data fromeach sensor is descriptive of waveform phenomena detected at the sensor.The method also includes determining an estimated spatial spectrum ofthe waveform phenomena based at least partially on the sensed data. Themethod further includes determining an estimated covariance matrix ofthe waveform phenomena based on the estimated spatial spectrum. Themethod includes determining adaptive beamforming weights using theestimated covariance matrix.

In another particular embodiment, a non-transitory computer-readablemedium includes instructions that, when executed by a processor, causethe processor to receive sensed data from sensors of a sensor array,where data from each sensor is descriptive of waveform phenomenadetected at the sensor. The instructions also cause the processor todetermine an estimated spatial spectrum of the waveform phenomena basedat least partially on the sensed data. The instructions further causethe processor to determine an estimated covariance matrix of thewaveform phenomena based on the estimated spatial spectrum. Theinstructions also cause the processor to determine adaptive beamformingweights using the estimated covariance matrix.

In another particular embodiment, a system includes a sensor arrayincluding multiple sensors and a signal processor coupled to the sensorarray. The signal processor is configured to receive sensed data fromthe sensors of the sensor array, where data from each sensor isdescriptive of waveform phenomena detected at the sensor. The signalprocessor is also configured to determine an estimated spatial spectrumof the waveform phenomena based at least partially on the sensed data.The signal processor is further configured to determine an estimatedcovariance matrix of the waveform phenomena based on the estimatedspatial spectrum. The signal processor is configured to determineadaptive beamforming weights using the estimated covariance matrix.

The features, functions, and advantages that have been described can beachieved independently in various embodiments or may be combined in yetother embodiments, further details of which are disclosed with referenceto the following description and drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of a particular embodiment of a signalprocessing system, including an array processor and a sensor array;

FIG. 2 is a diagram of a particular embodiment of a plane wave impingingsensors of the sensor array of FIG. 1;

FIG. 3 is a flow diagram of a first particular embodiment of a method ofsignal processing;

FIG. 4 is a flow diagram of a particular embodiment of a method ofdetermining a covariance matrix using multi-taper spectral estimationwith harmonic analysis;

FIG. 5 is a flow diagram of a second particular embodiment of a methodof signal processing; and

FIG. 6 is a block diagram of an illustrative embodiment of a computingenvironment.

DETAILED DESCRIPTION

In a particular embodiment, a covariance matrix for a sensor arrayresponsive to waveform phenomena (also referred to herein as aspace-time process) may be determined based on locations of individualsensor element of the sensor array, directional response and noise ofthe individual sensor elements, and a spatial spectrum of the waveformphenomena. In this embodiment, the covariance matrix may have aparticular structure that can be exploited to improve adaptivebeamformer performance both in terms of snapshot deficient operation andin terms of robustness against correlated signal and interferenceenvironments. Such performance improvements may be particularlybeneficial for large aperture arrays employing large numbers of sensorelements when operating in non-stationary and multi-path environments.

In a particular embodiment, the covariance matrix is determined usingstructured covariance techniques referred to herein as a covariance fromspatial spectrum (CSS) method. In disclosed CSS methods, the spatialspectrum of the sensor array observing the waveform phenomena isestimated, and spatial spectrum to covariance transforms are employed toestimate the covariance matrix. Performance predictions indicate thatthe disclosed structured covariance techniques can provide near optimalperformance for passive signal detection and recovery with very fewsnapshots (e.g., with a single snapshot or with fewer snapshots thaninterferers, as described further below).

In a particular embodiment, multi-taper spectral estimation withharmonic analysis is used to estimate the spatial spectrum. Otherspectral estimation techniques, such as classical power spectralestimation techniques may be used in other embodiments. The spatialspectrum estimation techniques disclosed may be extended to supportarbitrary array geometry and non-ideal array manifold response. Thenon-ideal array manifold response may include random ornon-deterministic variations in gain, phase, and directionality in thesensors, deterministic positional errors (e.g., the bending experiencedby an underwater towed array), or both non-deterministic anddeterministic errors.

FIG. 1 is a block diagram of a particular embodiment of a signalprocessing system 100. The signal processing system 100 includes asensor array 110 including multiple sensors 111-114. The sensor array110 may be coupled to an array processor 102. The sensor array 100 mayhave an arbitrary geometry. For example, the sensors 111-114 may beuniformly or non-uniformly spaced, arranged linearly or non-linearly(e.g., curved), and so forth.

The array processor 102 may be configured to receive sensed data fromthe sensors 111-114 of the sensor array 110. Data from each sensor111-114 may be descriptive of waveform phenomena detected at the sensor111-114. The array processor 102 may be configured to perform structuredcovariance estimation from spatial spectra for adaptive beamforming. Forexample, the array processor 102 may determine an estimated spatialspectrum of the waveform phenomena based at least partially on thesensed data. The array processor 102 may determine an estimatedcovariance matrix of the waveform phenomena based on the estimatedspatial spectrum. The array processor 102 may also determine adaptivebeamforming weights using the estimated covariance matrix. The arrayprocessor 102 or one or more other system components 106 may beconfigured to apply the adaptive beamforming weights to the sensed datato reconstruct a time-domain signal from the waveform phenomena.

In a particular embodiment, the array processor 102 uses a structuredcovariance determination method, which incorporates some a prioriknowledge or constraints. For example, the array processor 102 mayimpose a priori knowledge or constraints, such as an assumption that aspace-time process being observed is stationary in time, space, or both(e.g., a source of the waveform phenomena is stationary), or geometricinformation or constraints associated with the sensors 111-114 of thesensor array 110. Such a priori knowledge or constraints may reduce anumber of unknown quantities to estimate or a size of an allowablesolution space. Thus, structured covariance determination methods mayconvergence to a final solution with little data or reduced data in somescenarios.

In FIG. 1, the sensors 111-114 of the sensor array 110 spatially samplean environment and provide time varying output data to the arrayprocessor 102, which operates on the output data. The sampledenvironment may include signal(s), interference, noise or a combinationthereof. Certain of the signals and interference may originate fromdiscrete point sources, such as a point source signal 124 and a pointsource interference 126. Other signals or interference may be spatiallyspread fields, such as a spatially spread noise field 120 and aspatially spread signal 122.

The array processor 102 performs one or more operations using the outputdata from the sensor array 110. For example, the array processor 102 maydetermine whether a signal is present or not for a particular rangeand/or direction of arrival (which may be referred to as a “detectionproblem”). In another example, the array processor 102 may use theoutput data to reconstruct one or more components of a time-domainwaveform (which may be referred to as a “beamforming problem”).Knowledge of temporal characteristics of signals in the sampledenvironment may be exploited when such characteristics are known, as maybe the case, for example, in communications or active radar or sonarprocessing. Alternately, the signal may be unknown. In a particularembodiment, the array processor 110 uses multi-taper spectral estimationin a manner that combines non-parametric spectral estimation andharmonic analysis.

Adaptive beamforming may be used to suppress noise and interferencecomponents detected in the sampled environment and to enhance responseto desired signals. The array processor 102 may determine a vector ofbeamforming weights (referred to herein as a weighting vector, w) foruse in adaptive beamforming. The weighting vector may be used to combinethe sensor array 110 outputs, at a particular time, x _(m), to produce ascalar output, y_(m)=w ^(H) x _(m). When the covariance matrix, R=E{xx^(H) _(m)}, for the sensor array outputs is known (or estimated), anoptimal or near-optimal set of beamforming weights, w _(opt), may bedetermined as w _(opt)=γR ⁻¹ s, where s is a steering vector and γ is anormalization factor determined via an optimization process. Thecovariance matrix R may not be known a priori, and may therefore beestimated. When the noise and interference are Gaussian, a samplecovariance matrix, {circumflex over (R)}_(SCM), may be used as anunconstrained maximum likelihood (ML) estimate of R. The samplecovariance matrix may be determine as:

$\begin{matrix}{{\hat{\underset{\_}{R}}}_{SCM} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{\underset{\_}{x}}_{m}{\underset{\_}{x}}_{m}^{H}}}}} & (1)\end{matrix}$

Eqn. (1) is valid over a duration that the sampled environment can beconsidered stationary. Nonstationarity in the sampled environment mayarise in both sonar and radar applications and may be a limiting factorin the ability to estimate R from the output data of the sensor array110. In the case of sonar, a dominant noise mechanism in clutteredenvironments may be shipping. A combination of continuous motion ofsound sources relative to the sensor array 110 at short ranges andcorresponding small resolution cells of the sensor array 110 mayestablish a time scale over which a short-term stationarity assumptionis valid. In space-time adaptive processing (STAP) for radar, aresolution of range-Doppler cells may be susceptible to nonstationarityin the form of heterogeneous clutter, caused by ground clutter varyingover the surface of the earth. Nonstationary conditions may be addressedby using a limited number of snapshots of data from the sensor array110. A snapshot refers to a set of samples including no more than onesample from each sensor 111-114 of the sensor array 110. Snapshotdeficient operation refers to a condition in which a limited amount ofsnapshot data is available for processing. Snapshot deficient operationmay be a concern when the stationarity assumption applies only over alimit number of snapshots and the limit number of snapshots isinsufficient for the algorithm to converge. In addition to concernscaused by snapshot deficient operation due to nonstationarity,correlated signal and interference can also be problematic in adaptivebeamforming. Correlated interference may arise in multipath or smartjamming scenarios, resulting in signal cancellation effects.Applications that have no a priori knowledge of the inbound signal atthe desired direction of arrival, such as passive sonar, may be subjectto this effect.

In a particular embodiment, covariance from spatial spectrum (CSS)processing for uniform linear and regularly spaced array geometries maybe performed using power spectral estimation techniques (such aswindowed averaged periodogram techniques) and fast Fourier transform(FFT) processing. However, such techniques may estimate the covariancewith bias. Additionally, discrete or point source interference andspatial separation from the steering vector may degrade performance.Additional processing may be performed to detect, estimate, and subtractdiscrete components from such sources from available snapshot data priorto determining a final covariance matrix estimate. Such processingenables operation with negligible normalized signal to interference andnoise ratio (SINR) loss under a large range of conditions.

The signal processing system 100 may also include one or more othercomponents coupled to the array processor 102. For example, the one ormore other system components 106 may include devices that use thecovariance matrix, the beamforming weights, the sensed data, thetime-domain signal, or other data available to or generated by the arrayprocessor 102.

FIG. 2 is a diagram of a particular embodiment of plane wave impingingsensors 111-114 of the sensor array 110 of FIG. 1. In FIG. 2, the sensorarray 110 is an N element array with arbitrary sensor positions subjectto one or more incident plane waves 202. The location of an N^(th)sensor 114 may be specified in Cartesian coordinates as p _(n)=[x_(n),y_(n), z_(n)]^(T). Outputs of the sensors 111-114 resulting from theincident plane wave 202 may be described by an array manifold responsevector, v=((v_(n)(u)))_(n), where u is a unit directional vector 206that is opposite a direction of propagation 204 of the one or more planewaves 202.

Under a narrowband assumption, the effect of the plane wave 202 as itpropagates across the array 110 in space may be approximated by a phaseshift, and the array manifold response vector may be represented as:

$\underset{\_}{v} = \left( {\left( {\exp\left( {{- j}{\underset{\_}{k}}^{T}{\underset{\_}{p}}_{n}} \right)} \right)_{n} = \left( \left( {\exp\;\frac{j2\pi}{\lambda}{\underset{\_}{u}}^{T}{\underset{\_}{p}}_{n}} \right) \right)_{n}} \right.$

Each array element 111-114 may be sampled simultaneously at a time indexm to produce a snapshot, x _(m), where data elements of the snapshot areindividual sensor outputs. M snapshots may be produced for processing,one for each time index m. The snapshot x _(m) may include componentsthat derive from different sources. For example, the snapshot x _(m) mayinclude a contribution from point sources, which may be described as Va_(m), where V is a matrix of array manifold responses and a _(m) is avector of point source amplitudes. When only a single point source ispresent, the contribution from the point source may be described asva_(m), where v is the array manifold response and a_(m) is a scalarvalue of the individual point source amplitude. Both the vector of pointsource amplitudes, a _(m), and the scalar point source amplitude, a_(m),may be modeled as complex Gaussian. The snapshot x _(m) may also includea contribution from background or environmental noise, which may bedescribed as n _(b,m). The snapshot x _(m) may also include acontribution from internal noise of the sensors which may be describedas n _(w,m). Thus, the snapshot x _(m) may be a zero mean complexGaussian random vector, since it may be a linear combination ofuncorrelated, zero mean, complex Gaussian random vectors.

For the N element sensor array 110 of omni-directional sensors 111-114at locations p _(n), n=1, 2, . . . , observing waveform phenomena (e.g.,the plane wave 202) that can be described as a function f(t, p) withcorresponding covariance R_(f)(p), the covariance matrix for the arrayoutputs, ignoring sensor noise components, may include:R _(f)=((R _(f)(p _(r) −p _(c))))r,c

which can be expressed in matrix form using the frequency-wavenumberspectrum, P_(f)(k), and array manifold response vector, v(k).R _(f)=(2π)^(−C)∫^(. . .) ∫_(R) _(C) P _(f)( k ) v ( k ) v ^(H)( k )dk

where −C is the dimension of k. The waveform phenomena may also beinterpreted another way. A wavenumber restriction, |k|=2π/λ, impliesthat plane waves 202 observed by the sensor array 110 correspond toangles of arrival that may physically propagate to the sensor array 110.The waveform phenomena in this case may be modeled as the sum ofuncorrelated plane waves, that are distributed according to a normalizeddirectional distribution, G_(f) (ω, θ, φ), across all angles of arrivalto the sensor array 110, where ω is a radian frequency, θ is an anglefrom vertical in a spherical coordinate system, and φ is an azimuthangle in the spherical coordinate system.

FIG. 3 is a flow diagram of a first particular embodiment of a method ofsignal processing. The method may be performed by an array processor,such as the array processor 102 of FIG. 1, or by another signalprocessing system. At 302, sensed data may be receiving from sensors ofa sensor array. Data from each sensor may be descriptive of waveformphenomena detected at the sensor. For example, the waveform phenomenamay include a signal of interest and one or more interference signalsfrom one or more interferers. The sensed data may correspond to a singlesnapshot from the sensor array or more than one snapshot.

At 304, one or more estimated point source interference components maybe subtracted from the sensed data. For example, a harmonic analysistechnique may be used to subtract the point source interferencecomponents.

At 306, an estimated spatial spectrum of the waveform phenomena may bedetermined based at least partially on the sensed data. For example, theestimated spatial spectrum may include a frequency-wavenumber spectrumof the waveform phenomena or a direction-of-arrival spectrum of thewaveform phenomena. In a particular embodiment, the estimated spatialspectrum is determined by applying windowed, averaged periodogramanalysis to the sensed data, at 308. In another particular embodiment,the estimated spatial spectrum is determined by applying multi-taperspectral estimation to the sensed data, at 310.

At 312, an estimated covariance matrix of the waveform phenomena may bedetermined based on the estimated spatial spectrum. In a particularembodiment, a source of the waveform phenomena may be treated asstationary for purposes of estimating the covariance matrix. Theestimated covariance matrix may be determined based on informationregarding locations of the sensors of the sensor array. In a particularembodiment, the estimated covariance matrix converges with fewersnapshots than interferers.

In a particular embodiment, the estimated covariance matrix isdetermined as a combination of a first portion (e.g., a visible region)and a second portion (e.g., a virtual region). The first portion maycorrespond to one or more physically propagating waves of the waveformphenomena. The second portion may be associated with sensor noise. Inthis embodiment, the method may include determining an estimate of thesensor noise based at least partially on the second portion.

At 314, adaptive beamforming weights may be determined using theestimated covariance matrix. At 316, the adaptive beamforming weightsmay be applied to the sensed data to reconstruct a time-domain signalfrom the waveform phenomena. Thus, the covariance matrix may bedetermined in a manner that improves adaptive beamformer performance interms of snapshot deficient operation, in terms of robustness againstcorrelated signal and interference environments, and in terms of rate ofconvergence to a solution. Such performance improvements may beparticularly beneficial for large aperture arrays employing largenumbers of sensor elements when operating in non-stationary andmulti-path environments.

FIG. 4 is a flow diagram of a particular embodiment of a method ofdetermining a covariance matrix using multi-taper spectral estimationand with harmonic analysis. The method may be performed by an arrayprocessor, such as the array processor 102 of FIG. 1, or by anothersignal processing system. At 402, multi-taper spectral estimation may beapplied to sensed data from a sensor array (such as the sensor array 110of FIG. 1). Zero-padding may be used to increase sampling density of theestimated spatial spectrum. In a particular embodiment, a number oftapers used for the multi-taper spectral estimation is less than anumber of the sensors of the sensor array.

At 404, harmonic analysis is applied with the multi-taper spectralestimation to detect and subtract discrete components from the senseddata. Applying the multi-taper spectral estimation with the harmonicanalysis may include detecting discrete components of the waveformphenomena using outputs of the multi-taper spectral estimation, at 406.At 408, estimated discrete component parameters corresponding to eachdetected discrete component may be determined. At 410, an estimatedarray error vector corresponding to each detected discrete component maybe determined. At 412, a contribution of each detected discretecomponent may be subtracted from the sensed data to generate a residualcontinuous component. For example, the estimated array error vector maycorrespond to a linear minimum mean square error estimate of non-idealresponse error in the sensed data. At 414, an estimated continuouscomponent covariance matrix may be determined based on the residualcontinuous component.

In a particular embodiment, the estimated covariance matrix isdetermined as a numeric combination of the estimated continuouscomponent covariance matrix and an estimated discrete componentcovariance for each detected discrete component. In this embodiment, theestimated discrete component covariance for a particular detecteddiscrete component is determined based on the estimated array errorvector corresponding to the particular detected discrete component andbased on the estimated parameters of the particular detected discretecomponent.

FIG. 5 is a flow diagram of a second particular embodiment of a methodof signal processing. The method may be performed by an array processor,such as the array processor 102 of FIG. 1, or by another signalprocessing system. In particular, the method of FIG. 5 describes aprocess to determine a covariance matrix from spatial spectrum using amultitaper spectral estimation technique. At 502, an analysis region maybe defined. For example, defining the analysis region may includedetermining a number of tapers to be used for multitaper spectralestimation. At 504, a fast Fourier transform (FFT) process andzero-padding may be performed on sampled data from an array to formtapered outputs. Free parameter expansion may be used to facilitateforming the taper outputs. At 506, discrete components may be detectedwithin the sampled data. Harmonic analysis may be used to detect thediscrete components. The discrete components may be related to pointsource interference. At 508, discrete parameters, i.e., parameters ofthe discrete components (e.g., wavenumber, interference to noise ratio,etc.), may be estimated. At 510, the discrete components may be removedfrom the sample data, e.g., by subtracting contribution of the discretecomponents from the sample data. To illustrate, a coherent subtractionmethod or a null-projector method, as described further below, may beused to remove the discrete components. In a particular embodiment, themethod may return to 506 to detect additional discrete components. Thus,the detection, estimation of parameters and removal of discretecomponents may be iterated.

At 512, a continuous background spectrum may be determined or estimatedbased on the sample data with discrete components removed. At 514, acovariance matrix may be estimated based on the continuous backgroundspectrum.

Thus, a covariance matrix for a sensor array observing a stationaryspace-time process (also referred to herein as a waveform phenomena) maybe determined based on individual sensor element locations, directionalresponse and noise of those individual sensor elements, and a spatialspectrum of the space-time process. In a particular embodiment, thecovariance matrix may have a particular structure that can be exploited,improving adaptive beamformer performance both in terms of snapshotdeficient operation and robustness against correlated signal andinterference environments. The description below provides additionaldetail regarding determination of a covariance matrix and determinationof beamforming weights based on the covariance matrix.

Structured Covariance Techniques

Structured covariance methods incorporate some a priori knowledge of aproblem space to add additional constraints to the problem. Theseconstraints may include, for instance, a Toeplitz structure or blockToeplitz structure within the covariance matrix. The a priori knowledgemay be based on established models, such as the space-time processobserved being stationary, or the geometry of the array elements. Theseconstraints reduce the number of unknown quantities to estimate or thesize of the allowable solution space, and may convergence to a solutionwith very little data in some scenarios. Thus, they offer the potentialfor meaningful adaptive processor performance with lower snapshotrequirements than sample covariance methods and their reduced rankderivatives. Structured covariance algorithms may be applied to both thesnapshot deficient and correlated signal and interference problems.

Structured Covariance Based on Spatial Spectrum

An observed environment may be modeled as a spatially and temporallystationary space-time random process in many applications in arraysignal processing. This model of the space-time process can be describedusing a Cramér spectral representation theorem. This theorem provides adescription of a random process as a sum of uncorrelated plane waves,distributed across a visible region of an array as a function of azimuthand elevation or across the visible region in a frequency-wavenumberdomain. These domains are related to each other, and each properlydescribes the space-time process observed. A theory of second ordercharacterizations of space-time random processes shows that either ofthese descriptions (or representations) can be related to the space-timecovariance of the process at the output of the array. This is analogousto the relationship between power spectral density and auto-correlationin the case of time series analysis.

Because the space-time covariance can be related to either of theserepresentations, embodiments disclosed herein consider estimating thesespectral quantities first in order to compute the covariance at theoutput of the array. A process that uses estimated spectral quantitiesto estimate covariance is referred to herein as covariance from spatialspectrum (CSS).

A spectral representation model may include uncorrelated plane wavecomponents. When encountering correlated signal and interferencescenarios, processors may degrade in performance due to signalcancellation effects. Constraining the estimated covariance to astationary model may reduce the contribution of the correlation withinthe estimated covariance matrix, mitigating signal cancellation. Thestationary model may also impart a structure to the covariance matrix.Structured covariance matrix techniques may converge more rapidly thantraditional sample covariance matrix techniques and derivativetechniques. While there may be no simple closed form solution for amaximum likelihood (ML) estimate for a structured covariance for thegeneral problem of an unknown number of signals in non-white noise, anaturally intuitive interpretation of the process in theazimuth-elevation or frequency-wavenumber domains may be used to addressthe problem.

Spectral representation theorem in the context of angle of arrival tothe array considers the space-time process itself. Sensor noise may bemodeled as an uncorrelated noise component within each sensor and may bea part of the covariance observed by an adaptive processor. CSSprocessing for uniform linear and regularly spaced array geometries mayuse power spectral estimation techniques and Fast-Fourier Transform(FFT) processing.

The impact of discrete, or point source interference and its spatialseparation from the steering vector, s, may degrade performance.Additional steps may be used to detect, estimate, and subtract this typeof discrete component from available snapshot data prior to finalcovariance estimation. This allows operation with negligible normalizedSINR loss under a large range of conditions. Multi-taper spectralestimation technique, such as Thomson's multi-taper spectral estimationas described in Thomson, D. J. “Spectrum estimation and harmonicanalysis.” In Proceedings of the IEEE, Volume 70 (1982), 1055-1096,which may combine the convenience of non-parametric spectral estimationalgorithms and harmonic analysis techniques, work well in thisapplication. Performance assessment of the final technique was performedthrough simulation of various interference and noise environments, withexcellent normalized SINR loss performance obtained below an M=2Kthreshold for sample covariance type algorithms. Extensions to supportarbitrary array geometry are disclosed.

A correlated signal and interference environment generally cannot bemodeled as a stationary space-time process. The correlated componentsperturb the Toeplitz structure of the covariance matrix in the uniformlinear array (ULA) case. The expected value of the estimated covariancefor this problem, while constrained to an uncorrelated plane wavesmodel, contains a bias term related to the original correlated data.Impact of this bias term is negligible as array length increases.Simulation to compare sample covariance based techniques and CSStechniques, using an effective SINR metric, demonstrated successfulmitigation of signal cancellation effect in the CSS techniques. The CSStechniques mitigate the signal cancellation effect without loss ineffective aperture, and provide the same low snapshot supportrequirements as the stationary space-time process case. The CSStechniques provide an increase in effective sample size, similar toredundancy averaging, while maintaining a positive definite covariancematrix estimate, and differ from CMT in their ability to mitigate theeffects of correlated signal and interference in the data.

The CSS techniques are structured covariance techniques, which assume anideal array manifold response. A real-world array manifold response maybe non-ideal due to random variation in gain, phase, and directionalityin the sensors. Non-ideal response may also result from deterministicpositional errors, e.g., the bending experienced by an underwater towedarray. Performance impact for both random and deterministic arraymanifold response errors was simulated. In a particular embodiment, theCSS techniques are extended to account for non-ideal response.

Notation

The following notational conventions are used herein. Italicized lowercase letters (e.g., x) denote scalar variables, and x* is the complexconjugate of x. A vector of dimension P×1 is represented as a lowercaseletter with a single underline (e.g., k) and is complex or real-valuedas indicated. The transpose of a vector x is denoted as x ^(T), and x^(H) is the conjugate transpose of x. The matrix of dimension P× Q isdenoted with an underlined uppercase letter (e.g., X) and is complex orreal-valued as indicated. A compact notation is sometimes used to denotea matrix X=((x_(rc)))_(rc) and to indicate that the matrix X includeselements (rows, columns) which are individually denoted as x_(rc). Forexample, this notation may be used to denote values when the values area function of position indices, e.g., x=((cos [ωp]))_(p) is a vector ofP elements including the values [cos(ω), cos(2ω), cos(3ω)), •••]^(T).Similar notation may be used for vectors or higher dimensional matrices.The notation diag (x1, x2, •••) describes a diagonal matrix includingentries x₁, x₂, ••• on the main diagonal.

The expectation operator is denoted as E {x}. When describing randomvariables,

${x\overset{d}{\rightarrow}{{fx}(x)}},$indicates that the random variable x is distributed according to theprobability density function fx (x). Gaussian random variables aredenoted with specific notation for the appropriate density function. TheP-variant Gaussian random variable, x ε R^(P), with mean m=E {x} andcovariance R _(x)=E{(x−m)(x−m)^(T)} is represented using the notation

$\underset{\_}{x}\overset{d}{\rightarrow}{{{RN}_{P}\left( {\underset{\_}{m},{\underset{\_}{R}}_{x}} \right)}.}$Similarly, a P-variant complex Gaussian random variable, x ε C^(P), withmean m=E{x} and covariance R _(x)=E{(x−m)(x−m)^(H)} is represented usingthe notation

$\underset{\_}{x}\overset{d}{\rightarrow}{{{CN}_{P}\left( {\underset{\_}{m},{\underset{\_}{R}}_{x}} \right)}.}$It is assumed that all complex Gaussian random variables are circularlysymmetric of the type described by Goodman. The acronym i.i.d. standsfor independent and identically distributed.

Data Models and Assumptions

For ease of description, it is assumed herein that a sensed or observedwaveform phenomena (also referred to herein as an space-time process) isnarrowband and propagates in a homogeneous medium with velocity c andtemporal frequency f=c/λ. The space-time process is also assumed toinclude plane waves that are solutions to a homogeneous wave equation.This places a restriction on wavenumber k such that | k |=2π/λ, or saidanother way, plane waves correspond to physically propagating angles ofarrival to the sensor array. The space-time process may typically beconsidered to be spatially wide sense stationary (WSS) and zero mean.The spatially wide sense stationary property implies the spatialcovariance is a function of difference in position only. Additional, forWSS, plane wave components making up the process are uncorrelated. Thisrestriction is relaxed when dealing with the special condition ofcorrelated signal and interference, as this special condition violatesthe spatially stationary model. Snapshots are assumed independent over atime index, m. Several formulations of narrowband bandwidth and samplingperiod may be used to reasonably support this assumption.

General Model

The general model of an N element array with arbitrary sensor positionssubject to an incident plane wave is shown in FIG. 2. The nth elementlocation, specified in Cartesian coordinates, is given as p _(n)=[x_(n),y_(n), z_(n)]^(T). Sensor array element outputs resulting from incidentplane waves 202 are described by the array manifold response vector 216as v=((v_(n)(u)))_(n). For an array of omni-directional sensors (such asthe sensor array 110), under the narrowband assumption, the effect on aplane wave as it propagates across the array in space is approximated bya phase shift, and the array manifold response vector is v=((exp(−j k^(T) p _(n)))_(n)=((exp(j 2π/λ u ^(T) p _(n)))_(n). Each array elementis sampled simultaneously at time index m to produce a snapshot, x _(m)ε C^(N), where the elements of the snapshot vector are the individualsensor outputs, x _(m)=((x_(n)(m)))_(n). There are M snapshots availablefor processing.

In the output of the sensor array observing a space-time process, asnapshot model may include three terms

$\begin{matrix}{{\underset{\_}{x}}_{m} = {{\sum\limits_{k = 1}^{K}{{\underset{\_}{v}}_{k}{a_{k}(m)}}} + {\underset{\_}{n}}_{b,m} + {\underset{\_}{n}}_{w,m}}} & (2.1)\end{matrix}$

Each of the individual terms is described below. The space-time processmay include two types of sources. The first type corresponds to pointsources in the environment. These arrive at the sensor array 110 asdiscrete plane waves 202. K such sources may exist. Each source has agiven direction of arrival u _(k), with a corresponding array manifoldresponse vector, v _(k) as well as a source amplitude at each samplinginstant a_(k)(m). Combining the array manifold responses, V=[v ₁, v ₂,•••, v _(K)] and point source amplitudes a _(m)=((a_(k)(m)))_(k), Eqn.(2.1) can be written more succinctly asx _(m) =Va _(m) +n _(b,m) +n _(w,m)  (2.2)

The discrete source amplitudes, a _(m), are assumed to be

$\left. {i.i.d.\;{am}}\rightarrow{{{{CNK}\left( {0,{\underset{\_}{R}a}} \right)}{\underset{\_}{a}}_{m}}\overset{d}{\rightarrow}{{{CN}_{K}\left( {\underset{\_}{0},{\underset{\_}{R}}_{a}} \right)}.}} \right.$This is a standard model for passive sonar reception of far fielddiscrete sources. The point sources may be correlated with each other ornot. Specifying this attribute is a distinguishing feature betweenmodels. For a spatially stationary space-time process, the sources areuncorrelated by definition and R _(a)=diag (σ₁ ², σ₂ ², . . . σ_(K) ²).

The second component of the space-time process is the background orenvironmental noise. This noise is spread spatially and has a smooth,continuous distribution across some region of angle of arrival to thesensor array. The background noise at each snapshot,

${\underset{\_}{n}}_{\underset{\_}{b},m},{{{is}\mspace{14mu}{i.i.d.\;{\underset{\_}{n}}_{\underset{\_}{b},m}}}\overset{d}{\rightarrow}{{CN}_{N}\left( {\underset{\_}{0},{\underset{\_}{R}}_{b}} \right)}}$and is uncorrelated with the discrete components of the space-timeprocess.

Each sensor 111-114 also produces an internal noise, which isuncorrelated with internal noise of the other sensors and is independentacross snapshots. This noise is modeled as

${\underset{\_}{n}}_{\underset{\_}{w},m}\overset{d}{->}{{{CN}_{N}\left( {\underset{\_}{0},{\sigma_{w}^{2}\underset{\_}{I}}} \right)}.}$The sensor noise is uncorrelated with other components of the space-timeprocess.

Because it is a linear combination of uncorrelated, zero mean, complexGaussian random vectors, the snapshot x _(m) is a zero mean, complexGaussian random vector. It is distributed as

$\begin{matrix}{{{\underset{\_}{x}}_{m}\overset{d}{->}{{CN}_{N}\left( {\underset{\_}{0},{\underset{\_}{R}}_{x}} \right)}},{{\underset{\_}{R}}_{x} = {{{\underset{\_}{VR}}_{a}{\underset{\_}{V}}^{H}} + {\underset{\_}{R}}_{b} + {\sigma_{w}^{2}\underset{\_}{I}}}}} & (2.3)\end{matrix}$

Eqn. (2.3) is a model of general plane waves in non-white noise, and canrepresent line component, continuous, or mixed spectra type ofprocesses.

Space-Time Random Processes

Second (2^(nd)) order characterizations of WSS space-time randomprocesses may be used to develop linear array processors (such as thearray processor 102 of FIG. 1). A second order characterization mayprovide a complete specification a WSS space-time random process for aclass of Gaussian space-time processes. Under the assumptions statedabove, the WSS space-time random process can be described as a sum ofuncorrelated plane waves allocated according to a distribution acrossall angles of arrival to the sensor array.

Cramér spectral representation provides a description of a stochasticprocess in terms of an orthogonal process in a transform domain. Thespectral representation can be used to develop the more familiar secondorder characterizations for the process and their relationships, ratherthan by stated definition. It may become clear from the spectralrepresentation how a process may be characterized by one of thefollowing: i) a smooth, continuous power spectral density, ii) adiscrete (or line component) spectrum, or iii) a mixed spectrum,containing both smooth and line components. This formulation may beuseful in developing techniques to deal with mixed spectrum when theyare encountered, which may be beneficial for good adaptive beamformerperformance.

Second Order Characterization of Space-Time Random Processes

A zero mean, WSS complex space-time process, {f(t, p)}, is defined for−∞<t<∞, and over some dimensionality of Cartesian space,

^(C), typically

³ so that p=[p_(x), p_(y), p_(z)]^(T). The following relationshipsdefine second order central moments of the process.

Second Order Characterizations

Space-time covariance between two points Δp=p ₁−p ₂ and times τ=t₁−t₂ isdefined asR _(f)(τ,Δp )=E{f(t,p )f*(t−τ,p −Δp )}  (2.4)

The temporal frequency spectrum-spatial correlation function, alsoreferred to as the cross spectral density, is related to R_(f) (τ, Δp)by a Fourier transform of a time lag variable.

$\begin{matrix}{{S_{f}\left( {\omega,\underset{\_}{\Delta\; p}} \right)} = {\int_{- \infty}^{\infty}{{R_{f}\left( {\tau,\underset{\_}{\Delta\; p}} \right)}{\mathbb{e}}^{- {j\omega\tau}}{\mathbb{d}\tau}}}} & (2.5)\end{matrix}$

The frequency-wavenumber spectrum is related to S_(f) (ω, Δp) by aFourier transform of the Cartesian spatial coordinate Δp, with note ofthe reverse sign convention of the complex exponential. The wavenumbervector k has dimension C similar to Δp.

$\begin{matrix}{{P_{f}\left( {\omega,\underset{\_}{k}} \right)} = {\int{\ldots{\int_{\bullet^{c}}{{S_{f}\left( {\omega,\underset{\_}{\Delta\; p}} \right)}{\mathbb{e}}^{j{\underset{\_}{k}}^{T}\underset{\_}{\Delta\; p}}{\mathbb{d}\underset{\_}{\Delta\; p}}}}}}} & (2.6)\end{matrix}$

Because of the Fourier transform pair relationships, the followinginverse transforms may also be used.

$\begin{matrix}{\mspace{79mu}{{R_{f}\left( {\tau,\underset{\_}{\Delta\; p}} \right)} = {\left( {2\pi} \right)^{- 1}{\int_{- \infty}^{\infty}{{S_{f}\left( {\omega,\underset{\_}{\Delta\; p}} \right)}{\mathbb{e}}^{j\omega\tau}{\mathbb{d}\omega}}}}}} & (2.7) \\{\mspace{79mu}{{S_{f}\left( {\omega,\underset{\_}{\Delta\; p}} \right)} = {\left( {2\pi} \right)^{- C}{\int{\ldots{\int_{\bullet^{c}}^{\;}{{P_{f}\left( {\omega,\underset{\_}{k}} \right)}{\mathbb{e}}^{{- j}\;{\underset{\_}{k}}^{T}\underset{\_}{\Delta\; p}}{\mathbb{d}\underset{\_}{k}}}}}}}}} & (2.8) \\{{R_{f}\left( {\tau,\underset{\_}{\Delta\; p}} \right)} = {\left( {2\pi} \right)^{- 1}\left( {2\pi} \right)^{- C}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{j\omega\tau}{\mathbb{d}\omega}{\int{\ldots{\int_{\bullet^{c}}^{\;}{{P_{f}\left( {\omega,\underset{\_}{k}} \right)}{\mathbb{e}}^{{- j}\;{\underset{\_}{k}}^{T}\underset{\_}{\Delta\; p}}{\mathbb{d}\underset{\_}{k}}}}}}}}}} & (2.9)\end{matrix}$

These results may be specialized for the narrowband, independentsnapshot model. For a narrowband or monochromatic process at ω_(o),P_(f) (ω, k)=P_(f) (k)−2πδ(ω−ω_(o)). This simplifies the relationshipbetween the covariance and frequency-wavenumber spectrum.

$\begin{matrix}{{R_{f}\left( {\tau,\underset{\_}{\Delta\; p}} \right)} = {\left( {2\pi} \right)^{- C}{\mathbb{e}}^{{j\omega}_{0}\tau}{\int{\ldots{\int_{\bullet^{c}}^{\;}{{P_{f}\left( \underset{\_}{k} \right)}{\mathbb{e}}^{{- j}\;{\underset{\_}{k}}^{T}\underset{\_}{\Delta\; p}}{\mathbb{d}\underset{\_}{k}}}}}}}} & (2.10)\end{matrix}$

Further, under the assumption of independent snapshots, the covarianceis zero for τ≠0, further simplifying the relationship to

$\begin{matrix}{{{R_{f}\left( \underset{\_}{\Delta\; p} \right)} \equiv {R_{f}\left( {0,\underset{\_}{\Delta\; p}} \right)}} = {\left( {2\pi} \right)^{- C}{\int{\ldots{\int_{\bullet^{c}}^{\;}{{P_{f}\left( \underset{\_}{k} \right)}{\mathbb{e}}^{{- j}\;{\underset{\_}{k}}^{T}\underset{\_}{\Delta\; p}}{\mathbb{d}\underset{\_}{k}}}}}}}} & (2.11)\end{matrix}$

This shows that the narrowband, independent snapshot problem isprincipally a spatial problem, the temporal related aspects are notconsidered.

Matrix Representation

Consider an N element array of omni-directional sensors at locations p_(m) n=1, 2, . . . •••, observing a process represented as {f(t, p)}.The covariance matrix for the array outputs, ignoring sensor noisecomponents, is a matrix including elementsR _(f)=((R _(f)( p _(r) −p _(c))))_(r,c)  (2.12)

which can be expressed in matrix form succinctly using thefrequency-wavenumber spectrum and array manifold response vector, v(k).

$\begin{matrix}{{\underset{\_}{R}}_{f} = {\left( {2\pi} \right)^{- C}{\int{\ldots{\int_{\bullet^{c}}{{P_{f}\left( \underset{\_}{k} \right)}{\underset{\_}{v}\left( \underset{\_}{k} \right)}{{\underset{\_}{v}}^{H}\left( \underset{\_}{k} \right)}{\mathbb{d}\underset{\_}{k}}}}}}}} & (2.13)\end{matrix}$

Directional Distribution of Plane Waves

The space-time process may also be interpreted another way. Thewavenumber restriction, |k |=2π/λ, implies that plane waves observed bythe array correspond to angles of arrival that may physically propagateto the array, i.e., 0≦θ≦π and 0≦φ≦2π in spherical coordinates. Thestationary space-time process in this case may be modeled as the sum ofuncorrelated plane waves that are distributed according to a directionaldistribution, G_(f)(ω, θ, φ), across all angles of arrival to the array.G_(f) (ω, θ, φ) is similar to a probability density, and

$\begin{matrix}{{\frac{1}{4\pi}{\int_{0}^{\pi}{\int_{0}^{2\pi}{{G_{f}\left( {\omega,\theta,\phi} \right)}\sin\;\theta{\mathbb{d}\theta}{\mathbb{d}\phi}}}}} = 1} & (2.14)\end{matrix}$

The covariance may be related to the directional distribution. Considertwo sensors in the field at locations p+, p ₂. Each sensor has arespective frequency and directional response, H_(i)(ω, θ, φ), i=1, 2.The difference in position, Δp=p ₁−p ₂, can be described in Cartesian orspherical coordinates.

$\begin{matrix}{\underset{\_}{\Delta\; p} = {{{\underset{\_}{p\;}1} - {\underset{\_}{p}2}} = \left\{ \begin{matrix}\left\lbrack {{\Delta\; p_{x}},{\Delta\; p_{y}},{\Delta\; p_{z}}} \right\rbrack^{T} \\\left\lbrack {{s\;\sin\;\gamma\;\cos\;\zeta},{s\;\sin\;\gamma\;\sin\;\zeta},{s\;\cos\;\gamma}} \right\rbrack^{T}\end{matrix} \right.}} & (2.15)\end{matrix}$

The cross-spectral density between two sensors may be expressed inspherical coordinates as

$\begin{matrix}{{S_{f}\left( {\omega,s,\gamma,\zeta} \right)} = {\left( {4\pi} \right)^{- 1}{S_{f}(\omega)}{\int_{- \pi}^{\pi}{\int_{0}^{\pi}{{G_{f}\left( {\omega,\theta,\phi} \right)}{H_{1}\left( {\omega,\theta,\phi} \right)}{H_{2}\left( {\omega,\theta,\phi} \right)} \times {\exp\left( {j{\frac{\omega\; s}{c}\left\lbrack {{\sin\;{\theta sin}\;{{\gamma cos}\left( {\phi - \zeta} \right)}} + {\cos\;{\theta cos}\;\gamma}} \right\rbrack}} \right)}{\sin(\theta)}{\mathbb{d}\theta}{\mathbb{d}\phi}}}}}} & (2.16)\end{matrix}$

The function S_(f) (ω) in front of the integral translates the relativelevels specified by the directional distribution, G_(f) (·), to theabsolute levels observed at the array. The cross-spectral density andspace-time covariance are a Fourier transform pair.

$\begin{matrix}{{R_{f}\left( {\tau,s,\gamma,\zeta} \right)} = {\frac{1}{2\pi}{\int_{- \infty}^{\infty}{{S_{f}\left( {\omega,s,\gamma,\zeta} \right)}{\exp({j\omega\tau})}{\mathbb{d}\omega}}}}} & (2.17)\end{matrix}$

For a narrowband process, S_(f) (ω, s, γ, ζ)=S_(f) (s, γ,ζ)·2πδ(ω−ω_(o)). Computing the space-time covariance at τ=0, thecovariance and cross spectral density are seen to be the same.

$\begin{matrix}\begin{matrix}{{R_{f}\left( {0,s,\gamma,\zeta} \right)} = {\frac{1}{2\pi}{\int_{- \infty}^{\infty}{{S_{f}\left( {s,\gamma,\zeta} \right)}\bullet\; 2{{\pi\delta}\left( {\omega - \omega_{0}} \right)}{\exp({j\omega\tau})}{\mathbb{d}\omega}}}}} \\{= {S_{f}\left( {s,\gamma,\zeta} \right)}}\end{matrix} & (2.18)\end{matrix}$

Using Eqns. (2.16) and (2.18) with H_(i)(ω, θ, φ)=1, i=1, 2, thecovariance, R_(f)(s, γ, ζ) between two omnidirectional sensors in agiven noise field, G_(f) (ω, θ, φ), can be found. This is particularlyconvenient if the noise field has a known form that can be convenientlyrepresented in spherical or cylindrical harmonics. The covariance matrixfor an array of sensors is populated by values of this function, whereΔp=p _(r)−p _(c).R _(f)=((R _(f)(s,γ,ζ)))_(r,c)  (2.19)

Spectral Representation

Spectral representation theorem for stationary random processes, e.g.,Cramér spectral representation, provides a mechanism for describing theproperties of a WSS temporal or space-time random process in therespective transform domain, e.g., frequency or frequency-wavenumber.With this representation, the quantities of interest for a physicallygenerated process have intuitive interpretation and representation. Forexample, distribution of power may be related to the space-time randomprocess as a function of frequency (temporal or spatial). Definitionsand a description of properties of spectral representation are providedbelow. The discussion below also explains how spectral representationcan be used to develop the relationships between the more common secondorder moment characterizations of random processes. The spectralrepresentation and its properties for a continuous time random processare introduced, followed by a discussion of properties for a discretetime random process. Extensions for WSS space-time random processes, inlight of model assumptions, follows. In the following discussion, Ωdenotes a continuous time frequency variable (in radians/sec) and ωdenotes a discrete time frequency variable (in radians/sample).

Continuous Time

Let {f (t)} be a zero mean, complex-valued stationary random processdefined over time−∞<t<∞. Additionally, let {f (t)} be stochasticallycontinuous, a mild regularity condition that holds when theautocorrelation function, R (τ), is continuous at the origin. Then bythe Cramér spectral representation theorem, there exists an orthogonalprocess, {Z (Ω)}, such that for all t, a realization of the process,f(t), can be expressed as

$\begin{matrix}{{f(t)} = {\int_{- \infty}^{\infty}{{\mathbb{e}}^{{j\Omega}\; t}{\mathbb{d}{Z(\Omega)}}}}} & (2.20)\end{matrix}$

where equality is in the mean square sense and Eqn. (2.20) is aRiemann-Stieljes integral. The orthogonal process is zero meanE{dZ(Ω)}=0  (2.21)

and its corresponding covariance

$\begin{matrix}{{E\left\{ {{\mathbb{d}{Z\left( \Omega_{1} \right)}}{\mathbb{d}{Z^{*}\left( \Omega_{2} \right)}}} \right\}} = \left\{ \begin{matrix}0 & {\Omega_{1} \neq \Omega_{2}} \\{\mathbb{d}{S^{(I)}\left( \Omega_{1} \right)}} & {\Omega_{1} = \Omega_{2}}\end{matrix} \right.} & (2.22)\end{matrix}$

implies that disjoint frequencies are uncorrelated. The functionS^((I))(Ω), known as the integrated spectrum of {f (t)}, is a bounded,non-decreasing function. If S^((I))(Ω) is differentiable, then:

$\begin{matrix}{\frac{\mathbb{d}{S^{(I)}(\Omega)}}{\mathbb{d}\Omega} = {S(\Omega)}} & (2.23)\end{matrix}$

and the covariance then becomesE{dZ(Ω)dZ*(Ω)}=S(Ω)dΩ  (2.24)

The function S(Ω) is referred to as the power spectral density function.The integrated spectrum S^((I))(Ω) plays a role similar to a randomvariable cumulative distribution function, and many of the results fordistribution functions may be applied to it. In analog to the Lebesguedecomposition theorem for distribution functions, the integratedspectrum can be written as including three componentsS ^((I))(Ω)=S ₁ ^((I))(Ω)+S ₂ ^((I))(Ω)+S ₃ ^((I))  (2.25)

S^((I)) (Ω) is continuous, meaning its derivative exists for almost allΩ such that S^((I))(Ω)=∫_(−∞) ^(Ω)S(Ω′)dΩ′. A process with an integratedspectrum with only this term, S^((I))(Ω)=S₁ ^((I))(Ω), has a smoothbackground spectral component. White or colored noise processes (AR, MA,ARMA) are of this type, which is referred to herein as a process with apurely continuous spectrum.

S₂ ^((I)) (Ω) is a step function with jumps of size p_(l) at frequenciesΩ_(l), l=1, 2, . . . . A process with an integrated spectrum with onlythis term, S^((I))(Ω)=S₂ ^((I))(Ω), has a purely discrete spectrum orline spectrum. A harmonic random process,

${\left\{ {f(t)} \right\} = {\sum\limits_{l = 1}^{K}{A_{l}{\mathbb{e}}^{{j\Theta}_{l}}{\mathbb{e}}^{{j\Omega}_{l}t}}}},$where Θ_(l) is uniform [−π, π] and A_(l) is a real-valued randomvariable, has this type of spectrum.

S₃ ^((I))(Ω) is a continuous singular function. This type of pathologicfunction is not of practical use for spectral estimation, and may beassumed to be identically equal to 0 herein.

From this classification of spectra, and the decomposition described byEqn. (2.25), three types of stationary random processes may beencountered. First, those with a purely continuous spectrum. Secondthose with a purely discrete spectrum. And last, those which are acombination of both types, referred to herein as a mixed spectrumprocesses.

The autocorrelation function, R (τ), and power spectral density, S (Ω),of the process can be related through the spectral representation andthe properties of the orthogonal process. Starting with the definitionof the autocorrelation functionR(τ)=E{f(t)f*(t−τ)}  (2.26)

and replace f(t) with its spectral representation from Eqn. (2.20)

$\begin{matrix}{{R(\tau)} = {\int_{- \infty}^{\infty}{{\mathbb{e}}^{{j\Omega}_{l}t}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{{j\Omega}_{2}{({\tau - t})}}t}E\left\{ {{\mathbb{d}{Z\left( \Omega_{1} \right)}}{\mathbb{d}{Z^{*}\left( \Omega_{2} \right)}}} \right\}}}}}} & (2.27)\end{matrix}$

Using the uncorrelated property of the covariance at disjointfrequencies from Eqn. (2.22), Eqn. (2.27) simplifies to

$\begin{matrix}{{R(\tau)} = {{\int_{- \infty}^{\infty}{{\mathbb{e}}^{j\Omega\tau}{\mathbb{d}{S^{(I)}(\Omega)}}}} = {\int_{- \infty}^{\infty}{{S(\Omega)}{\mathbb{e}}^{j\Omega\tau}{\mathbb{d}\Omega}}}}} & (2.28)\end{matrix}$

where the result of the last integral assumes S^((I))(Ω) isdifferentiable. Because R (τ) is an integrable, deterministic function,Eqn. (2.28) indicates that R (τ) and S (Ω) form a Fourier transformpair, such that

$\begin{matrix}{{S(\Omega)} = {\frac{1}{2\pi}{\int_{- \infty}^{\infty}{{R(\tau)}{\mathbb{e}}^{- {j\Omega\tau}}{\mathbb{d}\tau}}}}} & (2.29)\end{matrix}$

Discrete Time

The spectral representation also applies to discrete time randomprocesses, {f [n] }, with modifications explained below. Given astationary, discrete time random process, {f [n] }−∞<n<∞, a realizationthe process, f [n], has spectral representation using the orthogonalprocess, {Z (ω)}, given by

$\begin{matrix}{{f\lbrack n\rbrack} = {\int_{- \pi}^{\pi}{{\mathbb{e}}^{{j\omega}\; n}{\mathbb{d}{Z(\omega)}}}}} & (2.30)\end{matrix}$

where equality is in the mean square sense. The limits of integrationare restricted to ±π to reflect the unambiguous range of the discretetime frequency. The orthogonal process is zero mean, E {dZ (ω)}=0, withcovariance

$\begin{matrix}{{E\left\{ {{\mathbb{d}{Z\left( \omega_{1} \right)}}{\mathbb{d}{Z^{*}\left( \omega_{2} \right)}}} \right\}} = \left\{ \begin{matrix}0 & {\omega_{1} \neq \omega_{2}} \\{\mathbb{d}{S^{(I)}(\omega)}} & {\omega_{1} = \omega_{2}}\end{matrix} \right.} & (2.31)\end{matrix}$

where the bounded, non-decreasing function S^((I))(ω) is the integratedspectrum of {f[n]}. The covariance implies disjoint frequencies areuncorrelated. If S^((I))(ω) is differentiable, then:

$\begin{matrix}{\frac{\mathbb{d}{S^{(I)}(\omega)}}{\mathbb{d}\omega} = {S(\omega)}} & (2.32)\end{matrix}$

and the covariance of the increments becomesE{dZ(ω)dZ*(ω)}=S(ω)dω  (2.33)

The function S (ω) is referred to as the power spectral densityfunction. Following the same procedure as in the continuous time case,the autocorrelation, R [l]=E {x (n) x*(n−1)} and integrated spectrum arerelated via the spectral representation by

$\begin{matrix}{{R\lbrack l\rbrack} = {{\int_{- \pi}^{\pi}{{\mathbb{e}}^{{j\omega}\; l}{\mathbb{d}{S^{(I)}(\omega)}}}} = {\int_{- \pi}^{\pi}{{S(\omega)}{\mathbb{e}}^{{j\omega}\; l}{\mathbb{d}\omega}}}}} & (2.34)\end{matrix}$

where the final expression assumes S^((I)) (ω) is differentiable. Thepower spectral density and autocorrelation form a Fourier transform pairsuch that

$\begin{matrix}{{S(\omega)} = {\frac{1}{2\pi}{\int_{l = {- \infty}}^{\infty}{{R\lbrack l\rbrack}{\mathbb{e}}^{{- {j\omega}}\; l}}}}} & (2.35)\end{matrix}$

where equality is in the mean square sense.

Sampling and Aliasing

A sampled WSS continuous time random process {f_(c) (t)}, −∞<t<∞, withuniform sampling period T, produces a WSS discrete random process {f_(d)[n]}, −∞<n<∞. The subscripts c and d are used to reinforce thedistinction between the continuous and discrete time domains. As astationary discrete time random process, f_(d) [n]=fc (t₀+nT) has aspectral representation

$\begin{matrix}{{f_{d}\lbrack n\rbrack} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{\mathbb{e}}^{{j\omega}\; n}{\mathbb{d}{Z_{d}(\omega)}}}}}} & (2.36)\end{matrix}$

where the orthogonal process {Z_(d) (ω)} has the spectral representationproperties described above. The autocorrelation and power spectraldensity form a Fourier transform pair.

$\begin{matrix}{{{R_{d}\lbrack l\rbrack} = {\int_{- \pi}^{\pi}{{S_{d}(\omega)}{\mathbb{e}}^{{j\omega}\; l}{\mathbb{d}\omega}}}},{{{Sd}(\omega)} = {\frac{1}{2\pi}{\int_{l = {- \infty}}^{\infty}{{R_{d}\lbrack l\rbrack}{\mathbb{e}}^{{- {j\omega}}\; l}}}}}} & (2.37)\end{matrix}$

The discrete time process second order characterizations, R_(d) [l] andS_(d) (ω), may be related to the continuous time counterparts, R_(c)(τ)and S_(c) (Ω). The discrete time autocorrelation function samples itscontinuous time counterpart. This follows directly from its definition

$\begin{matrix}\begin{matrix}{{R_{d}\lbrack l\rbrack} = {E\left\{ {{x_{d}\lbrack n\rbrack}{x_{d}^{*}\left\lbrack {n - l} \right\rbrack}} \right\}}} \\{= {E\left\{ {{x_{c}\left( {t_{0} + {nT}} \right)}{x_{c}^{*}\left( {t_{0} + {n\; T} - {lT}} \right)}} \right\}}} \\{= {R_{c}({lT})}}\end{matrix} & (2.38)\end{matrix}$

Two approaches may be taken to determine the form of the discrete timeprocess power spectral density. The first considers both R_(c) (τ) andS_(c) (Ω) as deterministic functions, and utilizes the properties of theFourier transform and the continuous time sampling indicated in Eqn.(2.38) to arrive at

$\begin{matrix}{{{S_{d}(\omega)} = {\frac{1}{T}{\sum\limits_{\beta = {- \infty}}^{\infty}{S_{c}\left( \frac{\omega - {2{\pi\beta}}}{T} \right)}}}},{{\omega } \leq \pi}} & (2.39)\end{matrix}$

The second approach utilizes the spectral representation directly toarrive at the same result.

The relationship in Eqn. (2.39) can be used to determine parameters forthe sampling period, T, such that the resultant discrete spectrum is anaccurate representation of the underlying continuous spectrum, i.e., theNyquist sampling criteria. However, consider an alternate interpretationof Eqn. (2.39). Even if the Nyquist sampling criteria is not met, usingan estimate of the discrete spectrum, although it may be an aliasedversion of the true spectrum, it is possible to determine correct valuesfor the covariance at the sample points indicated by Eqn. (2.38). Thismay be insufficient to completely reconstruct the underlying spectrum,but it enables estimation of covariance matrix values for adaptiveprocessing.

In practice a finite number of samples are available for processing.This is the effect of a finite duration observation window and has adirect impact on the ability to estimate the discrete spectrum, whichmay be addressed via formulation of a multitaper spectral estimator, asexplained below.

These concepts may be applied directly to uniform linear arrayprocessing problems. For regular arrays, such as arrays based on amultiple of a fixed spacing (e.g., minimum redundancy arrays),additional complexity is introduced due to the impact of the equivalentwindowing function on the ability to estimate the discrete spectrum.Arbitrary array design may also have this issue. Additionally, arbitraryarray design may add potentially unusual spectral combinations due tothe non-uniform spacing compared to the uniform structure described byEqn. (2.39).

Continuous Space-Time

The spectral representation extends to the case of a spatially andtemporally stationary multidimensional random process, {f (t, p)}. Givena realization of the process, f (t, p), there exists an orthogonalprocess, {Z (ω, k)}, such that for all t, p

$\begin{matrix}{{f\left( {t,\underset{\_}{p}} \right)} = {\int_{- \infty}^{\infty}\mspace{14mu}{\ldots\mspace{14mu}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{j{({{\omega\; t} - {{\underset{\_}{k}}^{T}\underset{\_}{p}}})}}{\mathbb{d}{Z\left( {\omega,\underset{\_}{k}} \right)}}}}}}} & (2.40)\end{matrix}$

The orthogonal process, {Z(ω, k)}, is zero mean E {dZ(ω, k)}=0, anduncorrelated across disjoint frequency-wavenumber bands.

$\begin{matrix}{{E\left\{ {{\mathbb{d}{Z\left( {\omega_{1},{\underset{\_}{k}}_{2}} \right)}}{\mathbb{d}{Z^{*}\left( {\omega_{3},{\underset{\_}{k}}_{4}} \right)}}} \right\}} = \left\{ \begin{matrix}0 & {{\left( {\omega_{1},{\underset{\_}{k}}_{2}} \right)\bigcap\left( {\omega_{3},{\underset{\_}{k}}_{4}} \right)} = 0} \\{{P\left( {\omega,\underset{\_}{k}} \right)}{\mathbb{d}\omega}{\mathbb{d}\underset{\_}{k}}} & {{\omega_{1} = \omega_{2}},{{\underset{\_}{k}}_{2} = {\underset{\_}{k}}_{4}}}\end{matrix} \right.} & (2.41)\end{matrix}$

where it is assumed that the integrated frequency-wavenumber spectrum,P_(x) ^((I)) (ω, k), is differentiable.

The relationships between covariance, cross spectral density, andfrequency-wavenumber spectrum may be derived from the multidimensionalspectral representation in Eqn. (2.40).

Directional Distribution of Plane Waves

As an alternative to the frequency-wavenumber domain in Eqn. (2.40), onemay define the orthogonal process across all angles of arrival on asphere.

$\begin{matrix}{{\mathbb{d}{Z\left( {\omega_{0},\underset{\_}{p}} \right)}} = {\frac{1}{4\pi}{\int_{0}^{\pi}{\sin\;\theta{\mathbb{d}\theta}{\int_{0}^{2\pi}{{\mathbb{e}}^{{- j}\; k_{0}{{\underset{\_}{a}}_{r}^{T}{({\theta,\phi})}}\underset{\_}{p}}{\mathbb{d}{Z\left( {\omega_{0},\theta,\phi} \right)}}}}}}}} & (2.42)\end{matrix}$

where

$\begin{matrix}{{f\left( {t,\underset{\_}{p}} \right)} = {\int_{- \infty}^{\infty}{{\mathbb{e}}^{{j\omega}\; t}{\mathbb{d}{Z\left( {\omega_{0},\underset{\_}{p}} \right)}}}}} & (2.43)\end{matrix}$

where k₀=2π/λ, and a _(r) (θ, φ) is a unit vector in the radialdirection. The orthogonal process, {dZ(ω₀, θ, φ)}, defines an integratedspectrum, S (ω) G^((I)) (ω, θ, φ), where G (·) is used to be consistentwith the discussion above of the second order characterization ofspace-time random processes. The function S (ω) scales the relativelevels defined in G (·) to the absolute levels seen at the array.Assuming G^((I)) (ω, θ, φ) is differentiable

$\begin{matrix}{{E\left\{ {{\mathbb{d}{Z\left( {\omega_{1},\theta_{1},\phi_{1}} \right)}}{\mathbb{d}{Z^{*}\left( {\omega_{2},\theta_{2},\phi_{2}} \right)}}} \right\}} = \left\{ \begin{matrix}0 & {{\theta_{1} \neq \theta_{2}},{\phi_{1} \neq \phi_{2}},{\omega_{1} \neq \omega_{2}}} \\{{S(\omega)}{G\left( {\omega,\theta,\phi} \right)}\frac{2\pi}{\sin\;\theta}{\mathbb{d}\omega}} & {{\theta_{1} = \theta_{2}},{\phi_{1} = \phi_{2}},{\omega_{1} = \omega_{2}}}\end{matrix} \right.} & (2.44)\end{matrix}$

The cross spectral density is defined by

$\begin{matrix}{{E\left\{ {{{dZ}\left( {\omega_{0},\underset{\_}{p}} \right)}{{dZ}^{*}\left( {\omega_{0},{\underset{\_}{p} - \underset{\_}{\Delta\; p}}} \right)}} \right\}} = {{S_{f}\left( {\omega_{0},\underset{\_}{\Delta\; p}} \right)}\frac{d\;\omega}{2\pi}}} & (2.45)\end{matrix}$

Relating Eqn. (2.45) to Eqn. (2.42), and using Eqn. (2.44):

$\begin{matrix}{{S_{f}\left( {\omega_{0},\underset{\_}{\Delta\; p}} \right)} = {\frac{1}{4\pi}{S(\omega)}{\int_{- \pi}^{\pi}{\int_{0}^{\pi}{{G_{f}\left( {\omega,\theta,\phi} \right)}{\mathbb{e}}^{{- j}\; k_{0}{{\underset{\_}{a}}_{r}^{T}{({\theta,\phi})}}\underset{\_}{\Delta\; p}}{\sin(\theta)}{\mathbb{d}\theta}{\mathbb{d}\phi}}}}}} & (2.46)\end{matrix}$

One may compare Eqn. (2.46) to the earlier expression in Eqn. (2.16),which also included directional response of the individual sensorelements and expanded the a _(r) ^(T) Δp terms in spherical coordinates.By defining the directional distribution, G_(f) (ω, θ, φ), and requiringdisjoint regions in angle space to be uncorrelated Eqn. (2.44) thisshows how the spectral representation underlies the model of thestationary space-time process as the sum of uncorrelated plane wavesdistributed over all directions of arrival to the array.

Optimal Beamforming

The following discussion reviews optimal beamforming techniques givenobservation of the interference and noise environment, and the relatedproblem when the desired signal is also present in the data. Normalizedsignal to interference and noise ratio (SINR) loss metric is a measureof the decrease in output SINR of an implemented beamformer compared toan optimal processor. The normalized SINR loss metric may be used toassess the performance of a given adaptive beamforming algorithm. Thismetric and its application are described below.

Minimum Variance Distortionless Response (MVDR)

Given a series of snapshots, x _(m) ε C^(N) with E{x _(m) x _(m) ^(H)}=R_(x), an adaptive processor may use a weighted linear combination of thesensor outputs to produce a scalar output signaly _(m) =w ^(H) x _(m)  (2.47)

The processor should pass the desired direction of arrival, specified bysteering vector s, undistorted. This constraint may be expressed asw ^(H) s=1  (2.48)

Expected output power of the processor isE{|y _(m)|² }=E{|w ^(H) x _(m)|² }=w ^(H) R _(x) w   (2.49)

The design criteria is to minimize the expected output power, subject tothe distortionless constraint. This design criteria is the same asmaximizing the output SINR. The optimization problem is thenarg _(w) min w ^(H) R _(x) w s.t. w ^(H) s=1  (2.50)

Using the method of Lagrange multipliers, the constrained cost functionto be minimized is thenJ( w )= w ^(H) R _(x) w +λ( w ^(H) s−1)+λ*( s ^(H) w−1)  (2.51)

The cost function is quadratic in w. Taking the complex gradient withrespect to w, and setting to zero

$\begin{matrix}{\frac{\partial{J\left( \underset{\_}{w} \right)}}{\partial{\underset{\_}{w}}^{*}} = {{{{\underset{\_}{R}}_{x}\underset{\_}{w}} + {\lambda\;\underset{\_}{s}}} = \left. 0 \right|_{\underset{\_}{w} = {\underset{\_}{w}}_{opt}}}} & (2.52)\end{matrix}$

the optimal set of weights isw _(opt) =−λR _(x) ⁻¹ s   (2.53)

Using Eqn. (2.53) in Eqn. (2.48) to solve for the Lagrange multiplier,gives λ=−(s ^(H) R _(x) s)⁻¹, and combining the two produces the finalweight vector

$\begin{matrix}{{\underset{\_}{w}}_{opt} = \frac{{\underset{\_}{R}}_{x}^{- 1}\underset{\_}{s}}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{x}^{- 1}\underset{\_}{s}}} & (2.54)\end{matrix}$

When snapshot data used to estimate R _(x) contains only a noise andinterference environment, this processor is referred to as a minimumvariance distortionless response (MVDR). In the event the desired signalis also present in the snapshot data, the same solution for the weightvector results; however, such a solution is sometimes referred to as aminimum power distortionless response (MPDR) to indicate the differencein the observed data. In practice, the distinction makes a significantdifference in terms of the required snapshot support to achieve goodperformance.

Normalized SINR Loss

For cases involving uncorrelated signal and interference, a standardmetric for performance of an adaptive beamformer is degradation in theoutput signal to interferer and noise ratio (SINR) compared to thatobtainable with an optimal processor. The normalized SINR loss, ξ, isdefined as

$\begin{matrix}{\xi = \frac{{SINR}_{a}}{{SINR}_{o}}} & (2.55)\end{matrix}$

The subscript o designates true quantities or optimal values, while thesubscript a designates the actual or estimated values. For convenience,normalized SINR loss can be expressed on a dB scale, as ξ_(dB)=−10 log₁₀(ξ). In this way, ξ_(dB)=1 implies an output SINR that is 1 dB lowerthan obtainable by an optimal processor. For the optimal processor, SINRis computed as

$\begin{matrix}{{SINR}_{o} = \frac{{{{\underset{\_}{w}}_{o}^{H}\underset{\_}{s}}}^{2}}{{\underset{\_}{w}}_{o}^{H}{\underset{\_}{R}}_{o}{\underset{\_}{w}}_{o}}} & (2.56)\end{matrix}$

while for an implemented processor, SINR may be computed as

$\begin{matrix}{{SINR}_{a} = \frac{{{{\underset{\_}{w}}_{a}^{H}\underset{\_}{s}}}^{2}}{{\underset{\_}{w}}_{a}^{H}{\underset{\_}{R}}_{o}{\underset{\_}{w}}_{a}}} & (2.57)\end{matrix}$

The general expression for ξ, not assuming a particular form for theweights, is therefore

$\begin{matrix}{\xi = {\frac{{{{\underset{\_}{w}}_{a}^{H}\underset{\_}{s}}}^{2}}{{\underset{\_}{w}}_{a}^{H}{\underset{\_}{R}}_{o}{\underset{\_}{w}}_{a}}\bullet\frac{{\underset{\_}{w}}_{o}^{H}{\underset{\_}{R}}_{o}{\underset{\_}{w}}_{o}}{{{{\underset{\_}{w}}_{o}^{H}\underset{\_}{s}}}^{2}}}} & (2.58)\end{matrix}$

For an adaptive beamformer with weights designed using the minimumvariance distortionless response (MVDR) criteria, using Eqn. (2.54) inEqn. (2.56), the SINR for the optimal processor isSINR_(o) =s ^(H) R ₀ ⁻¹ s   (2.59)

while using Eqn. (2.54) in Eqn. (2.57) yields the SINR for animplemented processor

$\begin{matrix}{{SINR}_{a} = \frac{\left( {{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{a}^{- 1}\underset{\_}{s}} \right)\left( {{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{a}^{- 1}\underset{\_}{s}} \right)}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{a}^{- 1}{\underset{\_}{R}}_{0}{\underset{\_}{R}}_{a}^{- 1}\underset{\_}{s}}} & (2.60)\end{matrix}$

Using Eqn. (2.59) and Eqn. (2.60) in Eqn. (2.55), the expression for ξbecomes

$\begin{matrix}{\xi = {\frac{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{a}^{- 1}\underset{\_}{s}}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{o}^{- 1}\underset{\_}{s}}\bullet\;\frac{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{a}^{- 1}\underset{\_}{s}}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{a}^{- 1}{\underset{\_}{R}}_{0}{\underset{\_}{R}}_{a}^{- 1}\underset{\_}{s}}}} & (2.61)\end{matrix}$

Eqn. (2.61) is a general expression for ξ when beamformer weights arefound via MVDR, but it does not give any insight into performance as itrelates to the quantities estimated for the underlying model. The matrixinverse operations also make it difficult to follow directly how modelparameters influence the performance, except under some simplifyingassumptions.

Spectral Estimation Techniques

In a particular embodiment, covariance for adaptive beamforming isestimated by first estimating the spatial or wavenumber spectrum. Thus,techniques that require an estimate of the covariance a priori, such asCapon's MVDR or MUSIC, are not usable for this purpose. Two maintechniques for spectral estimation based upon the data are described:windowed, averaged periodogram (a non-parametric technique) andmultitaper spectral estimation.

Classical (Nonparametric) Spectral Estimation

A windowed, averaged periodogram is a technique for spectral estimationfor uniform sampled data series. By applying a predetermined fixedwindow function or taper, w=((w[n]))_(n), to the data, the behavior ofthe spectral estimate can be controlled with regard to frequencyresolution and spectral leakage, e.g., sidelobe suppression. Thesequantities may be traded off against each other.

For simplicity, a single dimension time series or uniform linear arrayprocessing may be assumed. For either the time series or arrayprocessing problem the procedure is identical once the snapshots havebeen established. Time series processing may consider a contiguouscollection of N_(TOTAL) samples. This collection may then be subdividedinto M snapshots of N samples each. Within the time series, thesnapshots may be specified such that there is overlap of samples betweenadjacent snapshots. For the array processing application, each snapshotrepresents the simultaneous sampling of each of the N array elements. Ineither application, once obtained, each snapshot, x _(m), is multipliedelement-by-element with the windowing function, w [n].y _(m)=((x _(n) [n]w[n]))_(n) =x _(m)

w   (2.62)

where ⊙ is the element-by-element or Hadamard product. This windowedsnapshot data may then be Fourier transformed, e.g., using efficientfast Fourier transform (FFT) algorithms

$\begin{matrix}{{{Y_{m}\lbrack l\rbrack} = {\sum\limits_{n = 0}^{N - 1}{{y_{m}\lbrack n\rbrack}{\exp\left( {{- j}\;\frac{2\pi\;\ln}{N_{FFT}}} \right)}}}},{l = 0},1,\ldots\mspace{14mu},{L - 1}} & (2.63)\end{matrix}$

A value of N_(FFT) may be selected to more finely sample the underlyingspectrum, e.g., zero-padding, but the fundamental “resolution” of thetransform is constant based on the amount of available samples, N. Afinal estimated spectrum is the average, magnitude squared of theoutputs of Eqn. (2.63):

$\begin{matrix}{{\hat{P}\lbrack l\rbrack} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{Y_{m}\lbrack l\rbrack}}^{2}}}} & (2.64)\end{matrix}$

The frequency domain may be referred to as the transform domain of thetime series, although for array processing the wavenumber is the spatialfrequency, and once normalized by their respective sample period orsensor separation the two are equivalent. The averaged, modifiedperiodogram processing shown here using the FFT provides a fixedresolution and sidelobe (leakage) performance across the frequencydomain, based on the characteristics of the selected window. The arrayproblem encounters the non-linear mapping between physical angle ofarrival and wavenumber. If a fixed response across angle space isdesired, the window function becomes a function of angle of arrival, orw _(l)=((w_(l)[n]))_(n), such that Eqn. (2.63) becomes

$\begin{matrix}{{{Y_{m}\lbrack l\rbrack} = {\sum\limits_{n = 0}^{N - 1}{{w_{l}\lbrack n\rbrack}{x_{m}\lbrack n\rbrack}{\exp\left( {{- j}\frac{2\pi\;\ln}{N_{FFT}}} \right)}}}},{l = 0},1,\ldots\mspace{14mu},{N_{FFT} - 1}} & (2.65)\end{matrix}$

Multitaper Spectral Estimation (MTSE)

A multitaper algorithm formulates the spectral estimation problem as oneof estimating second order moments of a spectral representation of aprocess. The following discussion considers uniform sampled time series.Application of the method to uniform linear array processing inwavenumber space is immediate. For arbitrary array geometries oroperation in angle space, the concepts are the same though a specializedmultitaper design may be used and processing may be more computationallyintensive.

Spectral Estimation

In the following discussion the temporal “centering” term, e^(j(N-1)/2),has been omitted to simplify the discussion, and relationships are shownfor M>1 available snapshots. Given a stationary discrete random process,{x [n] }, −∞<n<∞, a realization of the process, x [n], may have aspectral representation

$\begin{matrix}{{x\lbrack n\rbrack} = {\int_{{- 1}/2}^{1/2}{{\mathbb{e}}^{j\; 2\pi\;{fn}}{\mathbb{d}{Z(f)}}}}} & (2.66)\end{matrix}$

where the covariance of the orthogonal increment process defines thepower spectral density.E{dZ(f)dZ*(f)}=S(f)df  (2.67)

The problem of spectral estimation is to estimate the covariance of thisprocess. However, dZ (f) is not observable directly from the available,limited samples x [n], 0≦n<N. While the impact of this data limitingoperation (or projection onto a finite number of samples) is obvious inthe time domain, its effect on the spectral representation of theprocess is less immediate. Taking the Fourier transform of the samples

$\begin{matrix}{{y(f)} = {\sum\limits_{n = 0}^{N - 1}{{x\lbrack n\rbrack}{\mathbb{e}}^{{- {j2}}\;\pi\;{fn}}}}} & (2.68)\end{matrix}$

and inserting Eqn. (2.66) into Eqn. (2.68) gives what may be referred toas the fundamental equation of spectral estimation

$\begin{matrix}{{y(f)} = {\int_{{- 1}/2}^{1/2}{\frac{\sin\; N\;\pi\left( {f - v} \right)}{\sin\;{\pi\left( {f - v} \right)}}{\mathbb{d}{Z(v)}}}}} & (2.69)\end{matrix}$

This result is a linear Fredholm integral equation of the first kind,and cannot be solved explicitly for dZ (v). This is in line with theinability to reconstruct the entire realization of the process, x [n],−∞<n<∞, from the limited sample observation. Eqn. (2.69) can be solvedapproximately, for a local region (f_(o)−W, f_(o)+W) using aneigenfunction expansion of the kernel

$\frac{\sin\; N\;{\pi(v)}}{\sin\;{\pi(v)}}$and a local least squares error criterion. The eigenfunction equation isgiven by

$\begin{matrix}{{{\lambda_{d}\left( {N,W} \right)}{Q_{d}\left( {N,W,f} \right)}} = {\int_{- W}^{W}{\frac{\sin\; N\;\pi(v)}{\sin\;{\pi(v)}}{Q_{d}\left( {N,W,v} \right)}}}} & (2.70)\end{matrix}$

where 0<W<½ is a design choice and N is the number of available samples.There are N solutions to Eqn. (2.70), indexed by the subscript d. Theeigenvalues, 0<λ_(d) (N, W)<1, give a measure of the concentration ofthe eigenfunction Q_(d) (N, W, f) within a desired region, [−W, W]. Forthis particular problem, the solutions to Eqn. (2.70) are known. TheQ_(d) (N, W, f) are the discrete prolate spheroidal wave functions(DPSWF), which are related to q_(d) (N, W, n), the discrete prolatespheroidal sequences (DPSS) by a Fourier transform

$\begin{matrix}{{Q_{d}\left( {N,W,f} \right)} = {ɛ_{d}{\int_{n = 0}^{N - 1}{{q_{d}\left( {N,W,n} \right)}{\mathbb{e}}^{{- j}\; 2\pi\;{fn}}}}}} & (2.71)\end{matrix}$

where ε_(d)=1 for d even, j for d odd. These sequences are also known asthe Slepian sequences. There are approximately 2 NW significanteigenvalues for these functions. Defining the Fourier transform of thewindowed samples, y_(m) ^((d))(f), as the d^(th) eigencoefficients ofthe data

$\begin{matrix}{{y_{m}^{(d)}(f)} = {\sum\limits_{n = 0}^{N - 1}{\frac{1}{ɛ_{d}}{q_{d}\left( {N,W,n} \right)}{x_{m}\lbrack n\rbrack}{\mathbb{e}}^{{- j}\; 2\pi\;{fn}}}}} & (2.72)\end{matrix}$

The d^(th) eigenspectra, Ŝ_(d)(f) may be is computed by averaging themagnitude squared of y_(m) ^((d)) (f) over all snapshots

$\begin{matrix}{{{\hat{S}}_{d}(f)} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{y_{m}^{(d)}(f)}}^{2}}}} & (2.73)\end{matrix}$

Due to the orthogonality of Q_(d) (N, W, f) over the interval [−W, W],for locally near flat spectra the eigenspectra are approximatelyuncorrelated. By averaging them, the overall variance of the finalestimate is improved. Before considering a method for combining theeigenspectra, one might be interested to use all N availableeigenfunctions to improve the variance by increased averaging, but thisis not recommended.

Looking at Multiple Tapers

To develop a better understanding of the difference between using amultitaper technique and a single taper classical technique, a simpleexample may be considered. For this example, consider three approachesto estimating the spectrum using a windowed technique. Averagingmultiple, uncorrelated estimates improves estimation accuracy. This canbe accomplished using non-overlapping Hann windows. Alternatively, amultitaper design uses overlapping but orthogonal windows. Both achieveimprovement due to averaging uncorrelated estimates, but one wouldexpect the multitaper design to perform better overall as itincorporates more of the sample data in each estimate. In order toimprove the resolution using the Hann window one may increase the lengthof the window (at the expense of providing fewer uncorrelatedestimates). In the extreme, a single Hann window has better resolutionthan the multitaper but achieves no improvement due to averaging.Alternate formulations of the Hann based approach may be used, such as50% overlap, etc., but in general the multitaper has better performancein terms of frequency resolution and overall improvement in variance dueto averaging.

There is a limit to the number of tapers that may be appliedmeaningfully, based on N and the W selected. As a rule of thumb, thereare 2 NW significant eigenvalues (sometimes more conservativelyestimated as 2 NW−1), indicating 2 NW tapers are highly concentrated inthe region [−W, W]. For example, a first few eigenfunctions may have amainlobe concentrated largely within [−W, W], but for higher numberedwindows the main lobe may be mostly outside the desired region, and thesidelobe level may increase substantially. This implies is that powerestimates based on these higher numbered windows may be heavilyinfluenced by frequency content outside the area of interest. Thiseffect may be referred to as a broad band bias, which may beundesirable, in particular for high dynamic range non-flat spectra.Limiting the number of employed windows such that D≦2 NW provides somerobustness against broad band bias automatically. Improvement may alsobe achieved by appropriately combining of the individual eigenspectra asdiscussed further below.

Combining Eigenspectra

With N given, and W and D specified, a multitaper method may computeindividual eigenspectra, Ŝ_(d) (f), using Eqn. (2.73). For an assumedflat spectrum, there is a fixed optimal weighting scheme for combiningthe individual Ŝ_(d) (f), however, this is of limited use. If one had apriori knowledge that the spectrum was white, altogether differentestimations techniques could be applied. For non-flat spectrum, adaptiveweighting schemes may reduce the contribution of broadband bias whilemaintaining estimation accuracy. An iterative method may be used todetermine an eigenspectra weighting function, h_(d) (f), for non-whitenoise. Begin with an initial estimate of the spectrum, Ŝ_(d) (f), usinga flat spectrum fixed weighting

$\begin{matrix}{{{\hat{S}}_{f}^{(0)}(f)} = \frac{\sum\limits_{d = 1}^{D}{\lambda_{d}{{\hat{S}}_{d}(f)}}}{\sum\limits_{d - 1}^{D}\lambda_{d}}} & (2.74)\end{matrix}$

where the superscript indexing is introduced to indicate the appropriateiteration. Estimating the variance of the process as

σ² = ∫_(−1/2)^(1/2)Ŝ_(f)⁽⁰⁾(f)𝕕f,the following iterative procedure may be performed.

$\begin{matrix}{{h_{d}^{(n)}(f)} = \frac{{\hat{S}}_{f}^{({n - 1})}(f)}{{\lambda_{d}{{\hat{S}}_{f}^{({n - 1})}(f)}} - {\left( {1 - \lambda_{d}} \right)\sigma^{2}}}} & (2.75) \\{{{\hat{S}}_{f}^{(n)}(f)} = \frac{\sum\limits_{d = 1}^{D}{{{h_{d}(f)}}^{2}{{\hat{S}}_{d}(f)}}}{\sum\limits_{d = 1}^{D}{{h_{d}(f)}}^{2}}} & (2.76)\end{matrix}$

Typically, only 3 iterations of Eqn. (2.75) followed by Eqn. (2.76) areused for convergence.

Free Parameter Expansion

The formulation for multitaper spectral estimation develops an estimateof the power spectral density at f₀ by considering a region [f_(o)−W,f_(o)+W]. This estimate is valid for any f in this region. For thisreason, for a specific Ŝ_(f) (f_(o)′), there remains a question of whichf_(o) in the range f_(o)′−W≦f_(o)≦f_(o)′+W is most appropriate as allare valid. This choice may be referred to as a free parameter expansionof f₀. The final Ŝ_(f)(f_(o)′) should be a weighted average of theseestimates, typically over a range |f_(o)′−f_(o)|≦0.8 W. In practice, theeigencoefficients in Eqn. (2.72) may be computed using FFT techniquesfor efficiency. Additional eigencoefficients for free parameterexpansion may be generated by continuing to use the FFT with additionalzero-padding. As used herein, the scalar multiplier NFPE refers to avalue such that the full zero-padded FFT size is N_(FFT)=N×NFPE. Thus,NFPE indicates the number of sampled “between points” available for freeparameter expansion averaging.

Harmonic Analysis

As explained above, a discrete random process may have a mixed spectrum,such that it includes two independent processes—one with a continuousspectrum and one with a discrete spectrum. In terms of its spectralrepresentation, a realization of a mixed spectrum stationary discreterandom process {x [n] } may have spectral representation

$\begin{matrix}{{x\lbrack n\rbrack} = {\int_{{- 1}/2}^{1/2}{{\mathbb{e}}^{j\; 2\pi\;{fn}}\left\lbrack {{\mathbb{d}{Z_{1}(f)}} + {\mathbb{d}{Z_{2}(f)}}} \right\rbrack}}} & (2.77)\end{matrix}$

where dZ₁ (f) corresponds to the continuous spectrum process and hasincrements in continuum from −½≦f≦½, while dZ₂ (f) corresponds to theline spectrum process and has increments only at the discrete locationsof the frequencies in the harmonic process, f_(k) for k=1, 2, •••, K.The line components (impulses) in the spectrum, which are caused by aportion of the process being a harmonic random process, causedifficulties with both classical and multitaper techniques. This is aresult of the modulation property of the Fourier transform. Windowing ofthe data in the time domain results in convolution in the frequencydomain. If a line component has large SNR, the result is unintendedspectral leakage across frequency. A harmonic analysis approach may beapplied to deal with this phenomenon.

At each frequency f_(o), the multiple tapers define a region in thefrequency domain, [f_(o)−W, f_(o)+W]. The continuous spectrum portion ofthe process, {dZ₁ (f)}, is non-zero throughout the region[f_(o)−W,f_(o)+W]. For now, assume a single line component may exist inthis region. If a line component exists at f_(o), the increment process{dZ₂(f)} is only non-zero at f_(o) within [f_(o)−W, f_(o)+W]. Eachrealization of the process, dZ₂ (f), provides a complex valued constantat f_(o). An analysis of variance (ANOVA) test may be used to detect thepresence of the potential line component. For a single snapshot, thisdetection problem is termed the constant false alarm rate (CFAR) matchedsubspace detector. Within the region [f_(o)−W, f_(o)+W], a subspacelocated at f_(o) only may be defined. A vector, q _(o), may be definedthat includes the mean value of each of the tapers

$\begin{matrix}{{\underset{\_}{q}}_{o} = \left\lbrack {{\sum\limits_{n = 0}^{N - 1}{q_{1}(n)}},{\sum\limits_{n = 0}^{N - 1}{q_{2}(n)}},\ldots\mspace{14mu},{\sum\limits_{n = 0}^{N - 1}{q_{D}(n)}}} \right\rbrack^{T}} & (2.78)\end{matrix}$

and use this to form a projection matrix, P _(q), for the subspaceP _(q) =q _(o)( q _(o) ^(H) q _(o))⁻¹ q _(o) ^(H)  (2.79)

The null projector for the subspace, P _(q) ^(l) ,P _(q) ^(l) =I−P _(q)  (2.80)

defines “everywhere else” in the region [f_(o)−W, f_(o)+W]. Forming avector of the eigencoefficients, y _(m) (f_(o))y _(m)(f _(o))=[y _(m) ⁽¹⁾(f ₀),y _(m) ⁽²⁾(f ₀), . . . ,y _(m) ^((D))(f₀)]^(T)  (2.81)

a detection statistic may be formed by taking a ratio of the power ofthe eigencoefficients in the f_(o)*subspace to the power of theeigencoefficients outside that subspace. Formally, the detectionstatistic F (f_(o)) may be computed as

$\begin{matrix}{{F\left( f_{o} \right)} = {\frac{\sum\limits_{m = 1}^{M}{{{\underset{\_}{y}}_{m}^{H}\left( f_{o} \right)}{\underset{\_}{P}}_{q}{{\underset{\_}{y}}_{m}\left( f_{o} \right)}}}{\sum\limits_{m = 1}^{M}{{{\underset{\_}{y}}_{m}^{H}\left( f_{o} \right)}{\underset{\_}{P}}_{q}^{\bot}{{\underset{\_}{y}}_{m}\left( f_{o} \right)}}}\underset{H_{0}}{\overset{H_{1}}{\gtrless}}{\Upsilon\;{TH}}}} & (2.82)\end{matrix}$

and compared to an appropriate threshold, γTH.

The importance of the detection of the line components in the spectrumis that they may be identified, and after estimation of their unknownparameters (amplitude, frequency, and phase) subtracted from theoriginal data. The residual sample data may then be subject to thespectral estimation algorithm, now with line components removed. Thefinal spectral estimate may numerically “add” the line components to thecontinuous spectrum, with appropriate scaling for SNR and estimationaccuracy.

Structured Covariance Matrices Based on Frequency Wavenumber SpectralEstimation

The discussion below explains how the frequency-wavenumber spectrum canbe used as a basis for covariance matrix estimation for WSS space-timeprocesses. First, the sensitivity of a model based adaptive beamformerto errors in estimates of the model parameters is considered.Performance of these processors with respect to the estimation ofindividual model parameters, such as interferer angle-of-arrival orinterferer to noise ratio (INR), is investigated. Second, relationshipsbetween the model for the array processing problem and the desiredcovariance are considered. In addition to the contribution fromphysically propagating waves seen by the array, contributions due to thesensor noise component may be included to maintain robustness. Thus, amodel may be defined that includes both of these elements, and thatproduces the proper covariance matrix estimate after transformation tothe space-time domain. This leads to an approach for covarianceestimation from spatial spectrum (CSS) for uniform arrays, and guidessubsequent discussion of CSS for arbitrary arrays. Concentrating on theuniform linear array case, the classical power spectral estimationtechniques discussed above are applied to the problem. Expected valuesfor the resulting covariance are developed and used to predictnormalized SINR loss performance. These predictions indicate that someform of harmonic analysis to mitigate the effects of line components inthe spectrum may be used. This allows the processor to maintainoperation in the low normalized SINR loss region for a wide range ofconditions.

Sensitivity of a Model Based Beamformer

In the following discussion, performance of a model-based adaptivebeamformer is assessed relative to an optimal beamformer. Analysis isperformed to determine performance sensitivity to estimation errors forthe individual model components. This analysis gives insight into whatparameters within the model matter most in terms of impact to beamformerperformance without specifying the exact form of the processor. Boundsare developed indicating the estimation accuracy to obtain an acceptableperformance for an acceptable normalized SINR loss, ξ.

Single Plane Wave in Spatially White Noise

Consider the simple case of a single plane wave interferer in spatiallywhite noise. This provides a basic understanding of sensitivity to errorin estimation of the interferer INR or wavenumber, and develops a basisfor more complicated models. The model for the covariance matrix forthis problem isR _(o)=σ_(n) ² v _(o) v _(o) ^(H)+σ_(w) ² I   (3.1)

where v _(o)=((e^(−jk) ^(o) ^(T) ^(p) ^(n) ))_(n) is the array manifoldresponse for wavenumber k _(o). The unknown quantities are:

σ_(w) ², the variance of the uncorrelated sensor noise;

σ_(n) ², the variance of the plane wave interferer; and

k _(o), the wavenumber of the plane wave interferer.

Optimal beamformer weights for desired wavenumber, ks, specifyingsteering vector, s, may be determined using Eqn. (2.54). Given itssimple structure, there is a closed form expression for R _(o) ⁻¹. Usingthe matrix inversion lemma.

$\begin{matrix}{{\underset{\_}{R}}_{o}^{- 1} = {\sigma_{w}^{- 2}\left\lbrack {\underset{\_}{I} - {\frac{\frac{\sigma_{n}^{2}}{\sigma_{w}^{2}}}{1 + {\frac{\sigma_{n}^{2}}{\sigma_{w}^{2}}{\underset{\_}{v}}_{o}^{H}{\underset{\_}{v}}_{o}}}{\underset{\_}{v}}_{o}{\underset{\_}{v}}_{o}^{H}}} \right\rbrack}} & (3.2)\end{matrix}$

The projection matrix for the subspace spanned by a single vector, v, isdefined as P _(v)=v(v ^(H) v)⁻¹ v ^(H) which can be rearranged as v ^(H)vP _(v)=vv ^(H). Defining the quantity

$\begin{matrix}{{\beta_{o} = {{\underset{\_}{v}}_{o}^{H}{{\underset{\_}{v}}_{o}\left( {\frac{\sigma_{w}^{2}}{\sigma_{n}^{2}} + {{\underset{\_}{v}}_{o}^{H}{\underset{\_}{v}}_{o}}} \right)}^{- 1}}},} & \;\end{matrix}$the matrix inverse is represented asR _(o) ⁻¹=σ_(w) ⁻² [I−β _(o) P _(v) _(o) ]  (3.3)

Eqn. (3.3), which is similar to a projection matrix on the null space ofv _(o), is a weighted subtraction of the projection onto the range spaceof v_(o). β_(o) is a measure of the interferer to noise ratio, so theoperation of the optimal beamformer can be interpreted as projecting outthe interferer based on its relative strength against the spatiallywhite noise. The complete expression for the weight vector is

$\begin{matrix}{{\underset{\_}{w}}_{o} = \frac{\left\lbrack {\underset{\_}{I} - {\beta{\underset{\_}{P}}_{v_{o}}}} \right\rbrack\underset{\_}{s}}{{{\underset{\_}{s}}^{H}\left\lbrack {\underset{\_}{I} - {\beta\;{\underset{\_}{P}}_{v_{o}}}} \right\rbrack}\underset{\_}{s}}} & (3.4)\end{matrix}$

Note that the denominator is a scalar, and does not affect the shape ofthe beam pattern other than as a gain.

In a particular embodiment, a model based adaptive processor may knowthe form of the covariance, estimate the unknown parameters, and use theestimates to determine the adaptive weights using Eqn. (3.4). As usedbelow, the subscript o is used to denote a true or optimal quantity,while the subscript a denotes an estimated quantity. The estimatedquantities reflect the true value and the estimation error, σ_(w,a)²=σ_(w) ²+Δσ_(w) ², σ_(n,a) ²=σn2+Δσ_(n) ², k_(a)=k_(o)+Δk. Theestimated covariance matrix, R _(a), and the corresponding weightvector, w _(a), areR _(a)=σ_(n,a) ² v _(a) v _(a) ^(H)+σ_(w,a) ² I   (3.5)R _(a) ⁻¹=σ_(w,a) ⁻² [I−β _(a) P _(v) _(a) ]  (3.6)

$\begin{matrix}{{\underset{\_}{w}}_{a} = \frac{\left\lbrack {\underset{\_}{I} - {\beta_{a}{\underset{\_}{P}}_{v_{a}}}} \right\rbrack\underset{\_}{s}}{{{\underset{\_}{s}}^{H}\left\lbrack {\underset{\_}{I} - {\beta_{a}{\underset{\_}{P}}_{v_{a}}}} \right\rbrack}\underset{\_}{s}}} & (3.7)\end{matrix}$

The weight vector does not depend on the absolute values of σ_(w,a) ²and σ_(n,a) ², only their ratio. Thus, the following quantities may bedefined:

$\begin{matrix}{{{INR}_{o} \equiv \frac{\sigma_{n}^{2}}{\sigma_{w}^{2}}},{{{INR}_{a} \equiv \frac{\sigma_{n,a}^{2}}{\sigma_{w,a}^{2}}} = {{INR}_{o}{\bullet\Delta}\;{INR}}}} & (3.8)\end{matrix}$

where 0<ΔINR<1 indicates an underestimate, and ΔINR>1 indicates anoverestimate. The following discussion considers the sensitivity of ξ tounder/over-estimation of INR, and whether there is a range of Δk that isalso tolerant such that ξ is near unity.

Over/Under-Estimate of INR, ΔINR≠1

For this analysis let Δk=0, so that P _(v) _(a) =P _(v) _(o) , toconcentrate on the behavior w.r.t. ΔINR. β_(a) may be expressed as

$\begin{matrix}{\beta_{a} = \frac{N}{{\left( {\sigma_{n}^{2}/\sigma_{w}^{2}} \right)^{- 1}\left( {\Delta\;{INR}} \right)^{- 1}} + N}} & (3.9)\end{matrix}$

Four scenarios resulting in two different approximations for performanceare described in Table 3.1. An additional condition is considered in theupper right case of Table 3.1 that (INR)⁻¹(ΔINR)⁻¹≦1, or ΔINR≧(INR)⁻¹.This puts a limit on the amount of underestimation considered andeliminates the situation of gross underestimation of INR when theinterferer is strong. As seen in Table 3.1, the two approximations toconsider are: Case 1− β_(a)≈1 and Case 2− β_(a)≈0.

TABLE 3.1 INR estimation error cases, single plane wave interferer:overestimate INR underestimate INR ΔINR ≧ 1 ΔINR ≦ 1 INR > 1(INR)⁻¹(ΔINR)⁻¹ << 1, β_(a) ≈ 1 (INR)⁻¹(ΔINR)⁻¹ ≦ 1, β_(a) ≈ 1 INR < 1(INR)⁻¹(ΔINR)⁻¹ ≧ 1, β_(a) ≈ 1 (INR)⁻¹(ΔINR)⁻¹ >> 1, β_(a) ≈ 0

Case 1. Overestimate or Slightly Underestimated INR, βa≈1

The general expression for the SINR loss is approximately

$\begin{matrix}{\xi \approx \frac{N - {\left( {1/N} \right){{{\underset{\_}{s}}^{H}{\underset{\_}{v}}_{o}}}^{2}}}{N - {\left( {\beta_{o}/N} \right){{{\underset{\_}{s}}^{H}{\underset{\_}{v}}_{o}}}^{2}}}} & (3.10)\end{matrix}$

The performance depends on the number of elements, N, the INR, and thewavenumber separation. For a uniform linear array with element spacingd, with broadside as the desired steering vector, s, Eqn. (3.10) becomes

$\begin{matrix}{\xi \approx \frac{N - {\left( {1/N} \right){{\sin^{2}\left( {k_{o}{{dN}/2}} \right)}/{\sin^{2}\left( {k_{o}{d/2}} \right)}}}}{N - {\left( {{\sigma_{w}^{2}/\sigma_{n}^{2}} + N} \right)^{- 1}{{\sin^{2}\left( {k_{o}{{dN}/2}} \right)}/{\sin^{2}\left( {k_{o}{d/2}} \right)}}}}} & (3.11)\end{matrix}$

Defining the normalized frequency as the ratio of the operationalfrequency to the design frequency of the array,f_(norm)=f|f_(o),f≦f_(o), the expression can be written in terms of theangle of incidence to the array, θ. Note the normalized wavenumberk_(o)d=−πf_(norm) cos θ.

$\begin{matrix}{\xi \approx \frac{\begin{matrix}{N - {\left( {1/N} \right)\sin^{2}}} \\{\left( {\pi\; f_{norm}{\cos\lbrack\theta\rbrack}{N/2}} \right)/{\sin^{2}\left( {\pi\; f_{norm}{{\cos\lbrack\theta\rbrack}/2}} \right)}}\end{matrix}}{\begin{matrix}{N - {\left( {{\sigma_{w}^{2}/\sigma_{n}^{2}} + N} \right)^{- 1}\sin^{2}}} \\{\left( {\pi\; f_{norm}{\cos\lbrack\theta\rbrack}{N/2}} \right)/{\sin^{2}\left( {\pi\; f_{norm}{{\cos\lbrack\theta\rbrack}/2}} \right)}}\end{matrix}}} & (3.12)\end{matrix}$

As θ moves away from the mainlobe at broadside, this expression is ≈1.

Case 2. Largely Underestimated INR, β_(a)≈0

When INR is largely underestimated, the adaptive beamformer reduces tothe conventional beamformer. Normalized SINR loss is highly dependent onINR, as the processor takes no action to specifically suppress theinterferer. The expression for SINR loss for a uniform linear array inthis case is approximately

$\begin{matrix}{\xi \approx {\frac{N}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{o}^{- 1}\underset{\_}{s}}\bullet\;\frac{N}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{o}\underset{\_}{s}}}} & (3.13)\end{matrix}$

The normalized SINR loss for the conventional (non-adaptive) beamformeris:

$\begin{matrix}{\xi \approx {\frac{N}{N - {\left( {\beta_{o}/N} \right){{{\underset{\_}{s}}^{H}{\underset{\_}{v}}_{o}}}^{2}}}\bullet\;\frac{N}{N + {{\sigma_{n}^{2}/\sigma_{w}^{2}}{{{\underset{\_}{s}}^{H}{\underset{\_}{v}}_{o}}}^{2}}}}} & (3.14)\end{matrix}$

Large underestimation of the INR results in performance of thenon-adaptive beamformer, effectively failing to take any correctiveaction in the weight determination to null the interference. Thenormalized SINR loss for this case shows the strong dependence betweenSINR loss and INR. This is expected for a non-adaptive processor, thelarger the interferer the worse the performance.

Wavenumber Offset, Δk≠0

Now consider the impact of Δk≠0 on the performance. Assume that ΔINR=1,i.e., there is no error in estimating the INR, so that β_(a)=β_(o). Theexpression for the approximate SINR loss in this case is

$\begin{matrix}{\xi \approx \left\lbrack {{{\sigma_{n}^{2}/\sigma_{w}^{2}}\frac{{{{\underset{\_}{v}}_{o}^{H}\underset{\_}{P}\overset{\bot}{v_{a}}\underset{\_}{s}}}^{2}}{{{\underset{\_}{P}\overset{\bot}{v_{a}}\underset{\_}{s}}}^{2}}} + 1} \right\rbrack^{- 1}} & (3.15)\end{matrix}$

This expression may be used to determine a Δk that achieves anacceptable SINR loss, ξ_(OK). This can then be compared to the CramérRao bound to see how reasonable it is in terms of the imagined modelbased processor. This results in the general expression for the singleinterferer in uncorrelated white noise case.

$\begin{matrix}{{{{\underset{\_}{v}}_{o}^{H}\underset{\_}{P}\overset{\bot}{v_{a}}\underset{\_}{s}}}^{2} \leq {{{\underset{\_}{P}\overset{\bot}{v_{a}}\underset{\_}{s}}}^{2}\frac{1 - \xi_{OK}}{\xi_{OK}{\sigma_{n}^{2}/\sigma_{w}^{2}}}}} & (3.16)\end{matrix}$

Assumptions about the array geometry may be used to further simply thisexpression. For a uniform linear array (ULA), performance may beanalyzed in terms of Δk, with an assumed desired steering vectorcorresponding to broadside, k_(s)=0. Due to the choice of k_(s), theinterferer wavenumber k_(o) is the separation in wavenumber between thetwo. The following expressions can be derived through analysis ofnormalized SINR loss of a model based adaptive beamformer as a functionof estimation accuracies of the model parameters (e.g., wavenumber andinterferer to noise ratio (INR)):

$\begin{matrix}{{\Delta\;{kd}} \leq \left( {\frac{1}{{\overset{.}{D}}^{2}\left( {k_{o}d} \right)} \cdot \left\lbrack {N - {\frac{1}{N}\frac{\sin^{2}\left( {k_{o}{{dN}/2}} \right)}{\sin^{2}\left( {k_{o}{d/2}} \right)}}} \right\rbrack \cdot \frac{1 - \xi_{OK}}{\xi_{OK}{\sigma_{n}^{2}/\sigma_{w}^{2}}}} \right)^{1/2}} & (3.17) \\{{\overset{.}{D}({kd})} = {\frac{1}{2} \cdot \frac{\begin{matrix}{{N\;\sin\left( {{kd}/2} \right){\cos\left( {{kdN}/2} \right)}} -} \\{{\sin\left( {{kdN}/2} \right)}{\cos\left( {{kd}/2} \right)}}\end{matrix}}{\sin^{2}\left( {{kd}/2} \right)}}} & (3.18)\end{matrix}$

Eqn. (3.17) includes oscillations that indicate there are areas in kdspace that are more tolerant to the estimation error, Δk. Thesecorrespond to nulls in the conventional beam pattern that providesufficient attenuation against the interference. A smoother bound forthe expression that eliminates oscillations and neatly spans the lowervalues may be determined, resulting in a final smoothed expression:

$\begin{matrix}{{\Delta\;{kd}} \leq \left( {{\frac{2{{\sin\left( {k_{o}{d/2}} \right)}}}{N} \cdot \left\lbrack {N - {\frac{1}{N}\frac{\sin^{2}\left( {k_{o}{{dN}/2}} \right)}{\sin^{2}\left( {k_{o}{d/2}} \right)}}} \right\rbrack}\frac{1 - \xi_{OK}}{\xi_{OK}{\sigma_{n}^{2}/\sigma_{w}^{2}}}} \right)^{1/2}} & (3.19)\end{matrix}$

While Eqn. (3.17) produces oscillations, Eqn. (3.19) smoothly bounds thebottom values. To be conservative, Eqn. (3.19) may be used, althoughEqn. (3.19) may be too conservative as the separation between theinterferer and the desired steering vector approaches zero. Thiscorresponds to the interferer residing in the main lobe of thebeamformer.

Comparison to the Cramér Rao Bound

The spatial frequency/direction of arrival estimation accuracy specifiedby Eqn. (3.19) results in a prescribed amount of normalized SINR loss,ξ_(OK). This accuracy can be compared to the Cramér Rao (CR) bound forthe case of a single plane wave in noise. The CR bound for a givennumber of snapshots, M, is

$\begin{matrix}{{C_{CR}({kd})} = {{\frac{1}{M}\left\lbrack {\frac{1}{N\;{\sigma_{n}^{2}/\sigma_{w}^{2}}} + \frac{1}{\left( {N\;{\sigma_{n}^{2}/\sigma_{w}^{2}}} \right)^{2}}} \right\rbrack}\frac{6}{\left( {N^{2} - 1} \right)}}} & (3.20)\end{matrix}$

Model for Covariance

Estimating the covariance matrix may begin by considering a model thatincorporates the components that make up the covariance at the output ofan N element array for a narrowband, stationary space-time process, f(t, Δp). Using the spectral representation theorem, the space-timeprocess can be represented as a sum of uncorrelated plane wavesdistributed as function of angle of arrival to the array, G(θ, φ), orwavenumber, G(k). The corresponding wavenumber spectrum for the processis proportional to this normalized distribution, P_(f) (k)=αG(k), whereα accounts for scaling the relative levels defined by G(k) to theabsolute power level seen at the array.

As explained above, a stationary random process may include twouncorrelated components, one corresponding to a continuous spectrumprocess, and another corresponding to a discrete spectrum, i.e.,harmonic, process. Independent, white sensor noise adds a thirdcomponent to the array output covariance. The m^(th) snapshot containingthese components is

$\begin{matrix}\begin{matrix}{{\underset{\_}{x}}_{m} = {{\sum\limits_{k = 1}^{K}{{\underset{\_}{v}}_{k}{a_{k}(m)}}} + {\underset{\_}{n}}_{b,m} + {\underset{\_}{n}}_{w,m}}} \\{= {{\underset{\_}{v}}_{a_{m}} + {\underset{\_}{n}}_{b,m} + {\underset{\_}{n}}_{w,m}}}\end{matrix} & (3.21)\end{matrix}$

where v _(k) is the array manifold response vector at spatial frequencyk _(k). The covariance for this model includes three uncorrelated partsbased on these componentsR _(x) =VR _(a) V ^(H) +R _(b) +R _(w)  (3.22)

where

${{\underset{\_}{R}}_{a} = {E\left\{ {{\underset{\_}{a}}_{m}{\underset{\_}{a}}_{m}^{H}} \right\}}},{{\underset{\_}{R}}_{b} = {E\left\{ {{\underset{\_}{n}}_{b,m}{\underset{\_}{n}}_{b,m}^{H}} \right\}}},{{\underset{\_}{R}}_{w} = {{E\left\{ {{\underset{\_}{n}}_{w,m}{\underset{\_}{n}}_{w,m}^{H}} \right\}} = {\sigma_{w}^{2}\underset{\_}{I}}}},$and spatial stationarity requires the plane waves be uncorrelated, sothat R _(a)=diag(σ₁ ², σ₂ ², . . . , σ_(K) ²).

Thus,

$\begin{matrix}{{\underset{\_}{R}}_{x} = {{\sum\limits_{k = 1}^{K}{\sigma_{k}^{2}{\underset{\_}{v}}_{k}{\underset{\_}{v}}_{k}^{H}}} + {\underset{\_}{R}}_{b} + {\sigma_{w}^{2}\underset{\_}{I}}}} & (3.23)\end{matrix}$

Alternatively, grouping the terms for the space-time process, R _(f)=VR_(a) V ^(H)+R _(b), separately from the sensor noise component:R _(x) =R _(f) +R _(w)  (3.24)

The matrix R _(f) includes all terms that correspond to physicalpropagating waves and can be decomposed via eigendecomposition

$\begin{matrix}\begin{matrix}{{\underset{\_}{R}}_{f} = {{\underset{\_}{Q}}_{f}{\underset{\_}{\Lambda}}_{f}{\underset{\_}{Q}}_{f}^{H}}} \\{= {\sum\limits_{n = 0}^{N - 1}{\lambda_{f,n}{\underset{\_}{q}}_{f,n}{\underset{\_}{q}}_{f,n}^{H}}}}\end{matrix} & (3.25)\end{matrix}$

Depending on the particular form of the space-time process, f (t, Δp),and the array geometry, R _(f) may have rank N_(f)<N. In that event,some of the eigenvalues will be zero valued.

$\begin{matrix}{{\underset{\_}{R}}_{f} = {\sum\limits_{n = 0}^{N_{f} - 1}{\lambda_{f,n}{\underset{\_}{q}}_{f,n}{\underset{\_}{q}}_{f,n}^{H}}}} & (3.26)\end{matrix}$

Regardless of the rank of R _(f), the N eigenvectors Q _(f) form acomplete orthonormal set for the space C^(N×N). Using (3.25) in Eqn.(3.24), with Q _(f) Q _(f) ^(H)=I the matrix R _(x) can be expressed

$\begin{matrix}\begin{matrix}{{\underset{\_}{R}}_{x} = {{{\underset{\_}{Q}}_{f}{\underset{\_}{\Lambda}}_{f}{\underset{\_}{Q}}_{f}^{H}} + {\sigma_{w}^{2}{\underset{\_}{Q}}_{f}{\underset{\_}{Q}}_{f}^{H}}}} \\{= {{\sum\limits_{n = 0}^{N - 1}{\left( {\lambda_{f,n} + \sigma_{w}^{2}} \right){\underset{\_}{q}}_{f,n}{\underset{\_}{q}}_{f,n}^{H}}} + {\sum\limits_{n = N_{f}}^{N - 1}{\sigma_{w}^{2}{\underset{\_}{q}}_{f,n}{\underset{\_}{q}}_{f,n}^{H}}}}} \\{= {{\underset{\_}{Q}}_{f}{\underset{\_}{\Lambda}}_{x}{\underset{\_}{Q}}_{f}^{H}}}\end{matrix} & (3.27)\end{matrix}$

As evident in Eqn. (3.27), the white noise contribution from R _(w)guarantees that all the eigenvalues are non-zero, so that overallcovariance, R _(x),is full rank.

Estimating Visible Space Covariance, {circumflex over (R)} _(vs)

Consider the covariance associated with the space-time process only, R_(f). As explained above, the transform relationship between thefrequency-wavenumber spectrum, P_(f)(k), and the space-time covariance,R_(f), isR _(f)=(2π)^(−C)∫ . . . ∫_(vs) P _(f)( k ) v ( k ) v ^(H)( k )dk  (3.28)

where C in this expression is the dimension of the wavenumber used, C=1,2 or 3. In Eqn. (3.28), the range of integration is restricted to thevisible region for the array, corresponding to physical propagatingwaves arriving at the array with some azimuth, 0≦φ≦2π, and elevation,0≦θ≦π. For a given estimate of the visible region frequency-wavenumberspectrum, {circumflex over (P)}_(vs)(k), the corresponding covariance ofthe space-time process can be determined using Eqn. (3.28){circumflex over (R)} _(vs)=(2π)^(−C)∫ . . . ∫_(vs) {circumflex over(P)} _(vs)( k ) v ( k ) v ^(H)( k )dk   (3.29)

Spectral estimation techniques used to form the estimate {circumflexover (P)}_(vs) (k) may not be able to distinguish between thecontribution of the space-time process, f (t, Δp), and the sensor noisecomponent apparent within the visible space. Any basis projection orsteered beam power measurement technique will see both the content fromf (t, Δp) and the sensor noise.{circumflex over (P)} _(vs)( k )={circumflex over (P)} _(f)( k)+{circumflex over (σ)}_(w) ²  (3.30)

{circumflex over (R)} _(vs) from Eqn. (3.29) will contain the sum ofboth.{circumflex over (R)} _(vs)={circumflex over (R)} _(f)+{circumflex over(σ)}_(w) ²(2π)^(−C)∫ . . . ∫_(vs) v ( k ) v ^(H)( k )dk   (3.31)

Observe from Eqn. (3.31) that the contribution due to the sensor noisecomponent, when viewed only across the visible region, appears as anadditional isotropic noise in the environment.

Visible and Virtual Space

For certain array geometries, or when operating below the designfrequency, there may be a significant additional virtual space inaddition to the visible space available to the array. To define virtualspace, the subspace spanned by the columns of a matrix A may be denotedas span (A)≡“A”. As explained above, with a sensor noise componentpresent, R _(x) ε C^(N×N), and is full rank, thus “R _(x)”=“C^(N×N)”.The visible space for a particular geometry and operational frequencymay only occupy a subspace of C^(N×N). “R _(x)” then includes twosubspaces, one corresponding to a visible region, indicated withsubscript vs, and one corresponding to a virtual region, indicated withsubscript yr.< R _(x) >=<R _(vs) >+<R _(vr)>  (3.32)

Wavenumbers in the virtual space do not correspond to physicalpropagating waves. Use of Eqn. (3.29) directly as an estimate ofcovariance with failure to account for the sensor noise component withinthe virtual region subspace may lead to poor sidelobe behavior withinthe virtual region, and adaptive beamformers developed using thiscovariance alone may suffer an overall loss in SINR.

Regularly and Uniformly Spaced Array Geometry

A regularly spaced array geometry is described by an interelementspacing that is a multiple of a fixed quantity, d, in Cartesiancoordinates. Uniform arrays are a special case of this with theinterelement spacing simply the constant, d. The uniform linear array,in terms of spatial sampling, is analogous to uniform sampled timeseries and the results of stationary processes and Fourier transformpairs and properties may apply directly. The discussion belowconcentrates on the uniform linear array, with some discussion of theimplications of regularly spaced arrays. Extensions to higher dimensionprocessing are envisioned.

The narrowband space-time process, f (t, Δp), at frequency ω includesplane waves propagating in a homogeneous medium with velocity c. Thesewaves are solutions to the homogeneous wave equation, and areconstrained in wavenumber such that | k |=ω/c=2π/λ. This indicates thatthe frequency-wavenumber spectrum, P_(f)(ω, k), for this process existson the surface of a sphere in wavenumber space with radius |k|.

Consider an N element uniform linear array with design frequency ω_(o)(spacing d=λ_(o)/2) and sensor elements at locations on the z axis,p_(n)=(n−1) d for n=0, •••, N−1. With no ability to resolve spatialcomponents in the k_(x) or k_(y) direction, the frequency-wavenumberspectrum, P_(f)(ω, k), may be projected down onto the k_(z) axis. Afterprojection the spectrum maintains the strict bandlimiting to the range|k_(z)|≦2π/λ. Using sampling of random processes, as explained above:

$\begin{matrix}{{R_{f}\left( {\Delta\; p} \right)} = {\frac{1}{2\pi}{\int_{{- 2}{\pi/\lambda_{o}}}^{2{\pi/\lambda_{o}}}{{P_{f}\left( k_{z} \right)}{\mathbb{e}}^{{- {j\Delta}}\;{pk}_{z}}{\mathbb{d}k_{z}}}}}} & (3.33)\end{matrix}$

where Δp=ld for integer l. For operation below the array designfrequency, ω<ω_(o), the wavenumber spectrum P_(f) (k_(z)) is non-zeroonly over the range corresponding to the visible region of the array,|k_(z)|≦2π/λ, and the range of integration may be reduced to ∫_(−2π/λ)^(2π/λ). Now consider the uncorrelated sensor noise componentR _(w)(Δp)=σ_(w) ²δ(Δp)  (3.34)

Even though it does not correspond to a component of a physicallypropagating space-time process, a wavenumber-spectrum and covarianceFourier transform pair can be defined:

$\begin{matrix}{{P_{w}\left( k_{z} \right)} = {{\sum\limits_{m = {- \infty}}^{\infty}{{R_{w}\left( {\Delta\; p} \right)}{\mathbb{e}}^{{j\Delta}\;{pk}_{z}}}} = \sigma_{w}^{2}}} & (3.35) \\{{R_{w}\left( {\Delta\; p} \right)} = {\frac{1}{2\pi}{\int_{{- 2}{\pi/\lambda_{o}}}^{2{\pi/\lambda_{o}}}{{P_{w}\left( k_{z} \right)}{\mathbb{e}}^{- {{j\Delta}{pk}}_{z}}{\mathbb{d}k_{z}}}}}} & (3.36)\end{matrix}$

The difference between P_(f) (k_(z)) and P_(w) (k_(z)) is that P_(w)(k_(r)) is non-zero over the entire interval, |k_(z)|≦2π/λ_(o). Thecovariance for the output of the array of interest is the sum of the twoR _(x)(Δ_(p))=R _(f)(Δp)+R _(w)(Δ_(p))  (3.37)

The two wavenumber spectra, one for the space-time process and the otherfor the sensor noise component, can be added to produce a compositespectrumP _(x)(k _(z))=P _(f)(k _(z))+P _(w)(k _(z)),|k _(z)|≦2π/λ_(o)  (3.38)

so that the covariance is related via

$\begin{matrix}{{R_{x}\left( {\Delta\; p} \right)} = {\frac{1}{2\pi}{\int_{{- 2}{\pi/\lambda_{o}}}^{2{\pi/\lambda_{o}}}{\left\lbrack {{P_{f}\left( k_{z} \right)} + {P_{w}\left( k_{z} \right)}} \right\rbrack{\mathbb{e}}^{- {{j\Delta}{pk}}_{z}}{\mathbb{d}k_{z}}}}}} & (3.39)\end{matrix}$

For convenience, the expressions can be converted to normalizedwavenumber, ψ=−k_(z)d space. Defining the ψ spectrum as

$\begin{matrix}{{P_{\psi,\alpha}(\psi)} = \left. {\frac{1}{d}{P_{\alpha}(k)}} \right|_{{k = {{- \psi}/d}},{\alpha = x},f,w}} & (3.40)\end{matrix}$

the covariance sequence R_(l,x) (l) and ψ spectrum are related

$\begin{matrix}{{R_{l,x}(l)} = {{R_{x}({ld})} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{P_{\psi,x}(\psi)}{\mathbb{e}}^{{j\psi}\; l}{\mathbb{d}\psi}}}}}} & (3.41)\end{matrix}$

P_(ψ,x) (ψ) may be estimated over the entire range, |ψ|≦π, from thesnapshot data, x _(m), using classical power spectral estimationtechniques. By accepting a fixed resolution in ψ space (which isnon-uniform in angle space, θ) this can be done efficiently with FFTbased processing. Using the estimate {circumflex over (P)}_(ψ,x) (ψ) inEqn. (3.41)

$\begin{matrix}{{{\hat{R}}_{l,x}(l)} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{\hat{P}}_{\psi,x}(\psi)}{\mathbb{e}}^{{j\psi}\; l}{\mathbb{d}\psi}}}}} & (3.42)\end{matrix}$

As long as the limits of integration in Eqn. (3.42) are over the entirerange, ∫_(−π) ^(π), the covariance estimate will contain the appropriatecomponents for both the space-time process and the uncorrelated sensornoise.

The process itself may be bandlimited, and if the sensor noise componentis estimated in the transform domain, it should have the appropriateform such that it equates to {circumflex over (σ)}_(w) ² I within thefinal covariance matrix estimate.

Positive Definiteness

In order for the estimated covariance to have value for adaptivebeamforming, the covariance should be Hermitian, {circumflex over (R)}_(x)={circumflex over (R)} _(x) ^(H), and invertible. The Hermitianproperty and invertibility imply that the eigenvalues of {circumflexover (R)} _(x) are all real valued and greater than zero, or more simplythat the matrix {circumflex over (R)} _(x) is positive definite. Forarbitrary vector x, not equal to the null vector (x ^(H) x≠0), thematrix A is positive semi-definite ifx ^(H) Ax≧0  (3.43)

and is indicated notationally as A≧0. The matrix A is positive definiteifx ^(H) Ax>0  (3.44)

and is indicated notationally as A>0.

For simplicity of explanation, consider a regularly spaced linear array(though higher dimension regular arrays may follow similarly). From Eqn.(3.42), the covariance matrix estimate is based on the estimate of thewavenumber spectrum, given in normalized wavenumber space, ψ=−kd

$\begin{matrix}{{\hat{\underset{\_}{R}}}_{x} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{\hat{P}}_{\psi,x}(\psi)}{{\underset{\_}{v}}_{\psi}(\psi)}{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\mathbb{d}\psi}}}}} & (3.45)\end{matrix}$

The ψ-spectrum estimates, {circumflex over (P)}_(ψ,x) (ψ), may berestricted to be real-valued and greater than zero. This hasimplications in choice of algorithm but is a reasonable for a powerspectral estimator when the observed process has a white noisecomponent. With this restriction, {circumflex over (P)}_(ψ,x) (ψ) may beexpressed in the form{circumflex over (P)} _(ψ,x)(ψ)={circumflex over (P)}_(ψ,f)(ψ)+{circumflex over (σ)}_(w) ²  (3.46)

where {circumflex over (P)}_(ψ,f)(ψ) is the estimate of the space-timeprocess spectrum, {circumflex over (P)}_(ψ,x)(ψ)≧0, and {circumflex over(σ)}_(w) ² is the estimate of the sensor noise {circumflex over (σ)}_(w)²>0. Now

$\begin{matrix}\begin{matrix}{{{\underset{\_}{x}}^{H}{\underset{\_}{\hat{R}}}_{x}\underset{\_}{x}} = {{{\underset{\_}{x}}^{H}\left( {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{\left\lbrack {{{\hat{P}}_{\psi,f}(\psi)} + {\hat{\sigma}}_{w}^{2}} \right\rbrack{{\underset{\_}{v}}_{\psi}(\psi)}{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\mathbb{d}\psi}}}} \right)}\underset{\_}{x}}} \\{= {{{\underset{\_}{x}}^{H}\left( {{\underset{\_}{\hat{R}}}_{f} + {\underset{\_}{\hat{R}}}_{w}} \right)}\underset{\_}{x}}}\end{matrix} & (3.47)\end{matrix}$

The first term is

$\begin{matrix}{{{\underset{\_}{x}}^{H}{\underset{\_}{\hat{R}}}_{f}\underset{\_}{x}} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{\hat{P}}_{\psi,f}(\psi)}{{{\underset{\_}{x}}^{H}{{\underset{\_}{v}}_{\psi}(\psi)}}}^{2}{\mathbb{d}\psi}}}}} & (3.48)\end{matrix}$

Both quantities in the integral in Eqn. (3.48), {circumflex over(P)}_(ψ,f)(ψ) and |x^(H)v_(ψ)(ψ)|², are real-valued and greater than orequal to zero, therefore {circumflex over (R)} _(f)≧0. The second termmay be evaluated as

$\begin{matrix}{{{\underset{\_}{x}}^{H}{\hat{\underset{\_}{R}}}_{w}\underset{\_}{x}} = {{\hat{\sigma}}_{w}^{2}{{\underset{\_}{x}}^{H}\left\lbrack {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{\underset{\_}{v}}_{\psi}(\psi)}{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\mathbb{d}\psi}}}} \right\rbrack}\underset{\_}{x}}} & (3.49)\end{matrix}$

The integral in Eqn. (3.29) can be reduced to

$\begin{matrix}{{\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{\underset{\_}{v}}_{\psi}(\psi)}{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\mathbb{d}\psi}}}} = \underset{\_}{I}} & (3.50)\end{matrix}$

so thatx ^(H) {circumflex over (R)} _(w) x={circumflex over (σ)} _(w) ² x ^(H)x>0  (3.51)

and {circumflex over (R)} _(w) is positive definite. The sum of apositive semidefinite matrix and a positive definite matrix is positivedefinite, so {circumflex over (R)} _(x)>0 when {circumflex over(P)}_(ψ,f)(ψ) is real-valued and greater than zero.

Performance when Using Classical PSD Techniques

The analysis of performance of estimating covariance from spatialspectrum (CSS) may begin by considering estimates of the wavenumberspectrum found using classical spectral estimation techniques. For easeof explanation, a uniform linear array with spacing d=λ_(o)/2 isconsidered. For a given fixed window function (or taper),w=((w[n]))_(n), the windowed snapshot data isy _(m)=((x _(m) [n]w[n]))_(n) =x _(m)

w   (3.52)

The windowed data may be used to form an averaged windowed periodogramestimate of the spectrum. Writing the array manifold response vector, v_(k)(k_(z))=((e^(−jk) ^(z) ^(dn)))_(n) in ψ=−k_(z)d space, v_(ψ)(ψ)=((e^(jψn)))_(n), the estimate may be developed in two steps.First, compute the Fourier transform of the windowed data.

$\begin{matrix}{{{\underset{\_}{Y}}_{m}(\psi)} = {{\sum\limits_{n = 0}^{N - 1}{{y_{m}(n)}{\mathbb{e}}^{{- {j\psi}}\; n}}} = {{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\underset{\_}{y}}_{m}}}} & (3.53)\end{matrix}$

The final spectral estimated is the averaged, magnitude squared value ofthe Fourier transforms

$\begin{matrix}\begin{matrix}{{{\hat{P}}_{y}(\psi)} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{Y_{m}(\psi)}}^{2}}}} \\{= {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\underset{\_}{y}}_{m}{\underset{\_}{y}}_{m}^{H}{{\underset{\_}{v}}_{\psi}(\psi)}}}}}\end{matrix} & (3.54)\end{matrix}$

{circumflex over (P)}_(ψ)(ψ) is periodic in ψ with period 2π. The range|ψ|≦π may be referred to as the region of support. The visible region ofthe array, when operating at frequency f=c/λ, (f≦f_(o)), is restrictedto the range |ψ|≦π(λ_(o)/λ). As explained above, the remainder outsidethe visible region is referred to as the virtual region. The fixedwindow function, w, provides a fixed resolution, e.g., “bin width”,across ψ space. This allows {circumflex over (P)}_(ψ)(ψ) to be computedefficiently at several equal spaced locations throughout the supportedregion using FFT techniques. It is also possible to have the windowfunction vary as a function of ψ, represented as w _(ψ). In this way,one can design for fixed resolution in angle space. This results innon-uniform “bin width” in ψ space. Multi-taper spectral estimationtechniques lend themselves to this method of design. FFT techniques maynot be directly applicable, though, when the window function is notfixed so increased computational cost may be associated with theapproach.

An alternate perspective on the estimated {circumflex over (P)}_(ψ)(ψ)spectrum may be considered in terms of an auto-correlation sequence{circumflex over (ρ)}_(y)[n] defined by the windowed sensor outputs.Based on the Fourier transform property

$\begin{matrix}{{{Y_{m}(\psi)}}^{2} = {F\left( {\sum\limits_{\beta = 0}^{N - 1}{{y_{m}\lbrack\beta\rbrack}{y_{m}^{*}\left\lbrack {\beta - n} \right\rbrack}}} \right)}} & (3.55)\end{matrix}$

The sample autocorrelation per snapshot is

$\begin{matrix}{{{{\hat{\rho}}_{y,m}\lbrack n\rbrack} = {\sum\limits_{\beta = 0}^{N - 1}{{y_{m}\lbrack\beta\rbrack}{y_{m}^{*}\left\lbrack {\beta - n} \right\rbrack}}}},{0 \leq n < N},{{{\hat{\rho}}_{y,m}\left\lbrack {- n} \right\rbrack} = {{\hat{\rho}}_{y,m}^{*}\lbrack n\rbrack}}} & (3.56)\end{matrix}$

where the sequence y_(m)[β] has value in the range [0, N−1], and is zeroelsewhere. As a convention, {circumflex over (ρ)}[n] may be used torepresent a sample autocorrelation from the data itself, while reservingR[n] is to indicate an auto-correlation based on the ensemble, E {·}.The overall sample autocorrelation is the average over all snapshots

$\begin{matrix}{{{\hat{\rho}}_{y}\lbrack n\rbrack} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{\hat{\rho}}_{y,m}\lbrack n\rbrack}}}} & (3.57)\end{matrix}$

Using Eqns. (3.55), (3.56), and (3.57), in Eqn. (3.54):

$\begin{matrix}{{{\hat{P}}_{y}\lbrack\psi\rbrack} = {\sum\limits_{n = {- {({N - 1})}}}^{N - 1}{{{\hat{\rho}}_{y}\lbrack n\rbrack}{\mathbb{e}}^{{- {j\psi}}\; n}}}} & (3.58)\end{matrix}$

The estimated power spectrum and auto-correlation sequence are a Fouriertransform pair, with corresponding inverse transform relationship

$\begin{matrix}{{{\hat{\rho}}_{y}\lbrack n\rbrack} = {\left( {2\pi} \right)^{- 1}{\int_{- \pi}^{\pi}{{{\hat{P}}_{y}(\psi)}{\mathbb{e}}^{{j\psi}\; n}{\mathbb{d}\psi}}}}} & (3.59)\end{matrix}$

The covariance matrix for the array is formed from the {circumflex over(ρ)}_(y) [n] values in a Toeplitz structure{circumflex over (R)} _(y)=(({circumflex over (ρ)}_(y)[r−c]))_(r,c)  (3.60)

Expressed directly in matrix notation based on Eqn. (3.59) this can bealso be expressed as

$\begin{matrix}{{\hat{\underset{\_}{R}}}_{y} = {\left( {2\pi} \right)^{- 1}{\int_{- \pi}^{\pi}{{{\hat{P}}_{y}(\psi)}{{\underset{\_}{v}}_{\psi}(\psi)}{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\mathbb{d}\psi}}}}} & (3.61)\end{matrix}$

It is also useful when comparing related techniques to understand howthe formation of {circumflex over (R)} _(y) relates to the operationsused in the traditional sample covariance matrix. Eqn. (3.57) isequivalent in result to Eqn. (3.59), but operates directly on thesnapshot data in space-time domain. First, define a windowed samplecovariance matrix

$\begin{matrix}{{\underset{\_}{R}}_{\underset{\_}{w},{SCM}} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{\left( {{\underset{\_}{x}}_{m}\bullet\underset{\_}{w}} \right)\left( {{\underset{\_}{x}}_{m}\bullet\underset{\_}{w}} \right)^{H}}}}} & (3.62)\end{matrix}$

The classical sample covariance matrix,

${{\underset{\_}{R}}_{SCM} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{\underset{\_}{x}}_{m}{\underset{\_}{x}}_{m}^{H}}}}},$uses Eqn. (3.62) with w=1, the all one's vector. Showing the entries inthe matrix in Eqn. (3.62) explicitly:

$\begin{matrix}{{\underset{\_}{R}}_{\underset{\_}{w},{SCM}} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{\quad\quad}\left( \begin{matrix}\begin{matrix}{x_{m,{\lbrack 0\rbrack}}w_{\lbrack 0\rbrack}x_{m,{\lbrack 0\rbrack}}^{*}w_{\lbrack 0\rbrack}^{*}} \\{x_{m,{\lbrack 1\rbrack}}w_{\lbrack 1\rbrack}x_{m,{\lbrack 0\rbrack}}^{*}w_{\lbrack 0\rbrack}^{*}} \\{x_{m,{\lbrack 2\rbrack}}w_{\lbrack 2\rbrack}x_{m,{\lbrack 0\rbrack}}^{*}w_{\lbrack 0\rbrack}^{*}}\end{matrix} & \begin{matrix}{x_{m,{\lbrack 0\rbrack}}w_{\lbrack 0\rbrack}x_{m,{\lbrack 1\rbrack}}^{*}w_{\lbrack 1\rbrack}^{*}} \\{x_{m,{\lbrack 1\rbrack}}w_{\lbrack 1\rbrack}x_{m,{\lbrack 1\rbrack}}^{*}w_{\lbrack 1\rbrack}^{*}} \\{x_{m,{\lbrack 2\rbrack}}w_{\lbrack 2\rbrack}x_{m,{\lbrack 1\rbrack}}^{*}w_{\lbrack 1\rbrack}^{*}}\end{matrix} & \ldots & {x_{m,{\lbrack 0\rbrack}}w_{\lbrack 0\rbrack}x_{m,{\lbrack{N - 1}\rbrack}}^{*}w_{\lbrack{N - 1}\rbrack}^{*}} \\\vdots & \; & \ddots & \; \\{x_{m,{\lbrack{N - 1}\rbrack}}w_{\lbrack{N - 1}\rbrack}x_{m,{\lbrack 0\rbrack}}^{*}w_{\lbrack 0\rbrack}^{*}} & \; & \ldots & {x_{m,{\lbrack{N - 1}\rbrack}}w_{\lbrack{N - 1}\rbrack}x_{m,{\lbrack{N - 1}\rbrack}}^{*}w_{\lbrack{N - 1}\rbrack}^{*}}\end{matrix} \right)}}}} & (3.63)\end{matrix}$

Comparing Eqn. (3.56) to Eqn. (3.63) shows that {circumflex over(ρ)}_(y,m) [n] is the sum down each of the diagonals in the inner matrixin Eqn. (3.63), where the diagonals correspond to the numbered index nas main, sub, or super diagonal according to:

$\begin{matrix}\begin{pmatrix}0 & {- 1} & {- 2} & \begin{matrix}\ldots & {- \left( {N - 1} \right)}\end{matrix} \\1 & \bullet & \bullet & \bullet \\2 & \bullet & \; & \; \\\vdots & \bullet & \; & \; \\{N - 1} & \; & \; & \;\end{pmatrix} & (3.64)\end{matrix}$

By averaging over multiple snapshots, therefore, {circumflex over(ρ)}_(y) [n] is the sum down the diagonals of R _(w,SCM). The covariancematrix {circumflex over (R)} _(y) is then populated with entries from{circumflex over (ρ)}_(y) [n]. Going forward as a convention, thisoperation is referred to herein as a diagonal-sum-replace (DSR), with anotation indicating the operation as{circumflex over (R)} _(y)=DSR( R _(w,SCM))  (3.65)

The DSR operation acts in a linear fashion for addition of matrices A, Bε C^(N×N)DSR( A+B )=DSR( A )=DSR( B )  (3.66)

as well as in regards to the expectation operatorE{DSR( A )}=DSR(E{A})  (3.67)

Expected Value-Stationary Random Process

The expected value of the covariance, {circumflex over (R)} _(y), forthe WSS space-time process is considered below. From Eqn. (3.58), theexpected value of {circumflex over (P)}_(y) [ψ] is related to the sampleautocorrelation, {circumflex over (ρ)}_(y) [n],

$\begin{matrix}{{E\left\{ {{\hat{P}}_{y}\lbrack\psi\rbrack} \right\}} = {\sum\limits_{n = {- {({N - 1})}}}^{N - 1}{{\mathbb{e}}^{{- {j\psi}}\; n}E\left\{ {{\hat{\rho}}_{y}\lbrack n\rbrack} \right\}}}} & (3.68)\end{matrix}$

The expected value of E {{circumflex over (ρ)}_(y)}=E{{circumflex over(ρ)}_(y,m)}.

$\begin{matrix}{{E\left\{ {{\hat{\rho}}_{y,m}\lbrack n\rbrack} \right\}} = {{R_{x}\lbrack n\rbrack}{\sum\limits_{\beta = 0}^{N - 1}{{w\lbrack\beta\rbrack}{w^{*}\left\lbrack {\beta - n} \right\rbrack}}}}} & (3.69)\end{matrix}$

where R_(x)[n] is the ensemble covariance. The remaining summation termis the sample autocorrelation of the window,

${\rho_{w}\lbrack n\rbrack} = {\sum\limits_{m = 0}^{N - 1}{{w\lbrack m\rbrack}{{w^{*}\left\lbrack {m - n} \right\rbrack}.}}}$The final result for the expectation is thenE{{circumflex over (ρ)} _(y) [n]}=R _(x) [n]ρ _(w) [n]  (3.70)

From Eqn. (3.70), the expected value of the covariance matrix isE{{circumflex over (R)} _(y) }=R _(x)

R _(w)  (3.71)

where R _(w)=DSR(ww ^(H))=((ρ_(w)[r−c]))_(r,c). Looking at the result inthe frequency domain, using Eqn. (3.70) in Eqn. (3.68) results in:

$\begin{matrix}{{E\left\{ {{\hat{P}}_{y}(\psi)} \right\}} = {\sum\limits_{n = {- {({N - 1})}}}^{N - 1}{{R_{x}\lbrack n\rbrack}{\rho_{w}\lbrack n\rbrack}{\mathbb{e}}^{{- j}\;\psi\; n}}}} & (3.72)\end{matrix}$

This can be expressed in the ψ domain as the convolution of the powerpattern of the window, C_(w)(ψ)=|W(ψ)|²=F (ρ_(w)[n]), with theunderlying model spectrum, P_(x,ψ) (ψ)

$\begin{matrix}\begin{matrix}{{E\left\{ {{\hat{P}}_{y}(\psi)} \right\}} = {{C_{w}(\psi)}\#{P_{x,\psi}(\beta)}}} \\{= {\frac{1}{2\;\pi}{\int_{- \pi}^{\pi}{{C_{w}\left( {\psi - \beta} \right)}{P_{x,\psi}(\beta)}\ {\mathbb{d}\beta}}}}}\end{matrix} & (3.73)\end{matrix}$

Performance Based on Expected Value

The result for the expected value of the covariance of Eqn. (3.71) canbe used to assess performance of the algorithm. For a given windowfunction (a.k.a. taper), w, first determine the matrix R _(w)=DSR (ww^(H)). For each particular problem of interest, e.g., single plane wavein uncorrelated noise, form the known model ensemble covariance, R _(x).Using these in the normalized SINR loss expression, adaptive beamformerperformance using CSS can be analyzed.

$\begin{matrix}{\xi = {\frac{{{\underset{\_}{s}}^{H}\left( {{\underset{\_}{R}}_{x}\mspace{11mu}\bullet\mspace{11mu}{\underset{\_}{R}}_{w}} \right)}^{- 1}\underset{\_}{s}}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{x}^{- 1}\underset{\_}{s}} \cdot \frac{{{\underset{\_}{s}}^{H}\left( {{\underset{\_}{R}}_{x}\mspace{11mu}\bullet\mspace{11mu}{\underset{\_}{R}}_{w}} \right)}^{- 1}\underset{\_}{s}}{{{\underset{\_}{s}}^{H}\left( {{\underset{\_}{R}}_{x}\mspace{11mu}\bullet\mspace{11mu}{\underset{\_}{R}}_{w}} \right)}^{- 1}{{\underset{\_}{R}}_{x}\left( {{\underset{\_}{R}}_{x}\mspace{11mu}\bullet\mspace{11mu}{\underset{\_}{\; R}}_{w}} \right)}^{- 1}\underset{\_}{s}}}} & (3.74)\end{matrix}$

Because of the Hadamard product nature of the relationship severalthings become apparent:

-   -   1. As the window function becomes ideal, the corresponding        ρ_(w)[n]=1, and R _(w)=1, the Hadamard product identity (the all        1's matrix). The estimated covariance goes to the ideal,        {circumflex over (R)} _(y)→R _(x), and there is no performance        loss. This corresponds to an idealized power pattern        C_(w)(ψ)=2πδ(ψ). Such a power pattern requires an infinite        number of elements in the array. For practical arrays there will        be a bias to the resultant covariance matrix because of the        non-zero width of C_(w)(ψ), indicating that, in general, this        approach will not converge to the optimal solution. Of        particular interest is the impact of the bias, and whether there        are conditions where the impact to performance is negligible.    -   2. It is difficult to simplify analysis because of the        element-by-element nature of the operation. There is no easy        expression for the inverse of a Hadamard product, and other        properties that are available are not helpful in the analysis of        normalized SINR loss. Thus, defined scenarios may be considered        and simulated over the range of parameters, such as        (ψ_(s)−ψ_(o)) and σ_(n) ²/σ_(w) ² in the single plane wave        interferer case.    -   3. The particular form of the classical window functions used        for spectral estimation result in R _(w) that is a normalized        diagonally homogeneous correlation matrix, meaning it has a        constant value on its main diagonal. From the eigenvalue        majorization theorem, the eigenvalue spread of R _(x) ⊙R _(w) is        less than or equal to the eigenvalue spread of R _(x). One can        also arrive at the same conclusion by considering that: 1) the        eigenvalues are bounded by the min/max of the power spectrum;        and 2) the technique used here effectively convolves the        underlying power spectrum with a smoothing window (the power        pattern), which will decrease the resultant min/max of the        smoothed spectrum compared to the original. This property may        lower the condition number of the matrix, cond

${\left( \underset{\_}{A} \right) = \frac{\lambda_{\max}\left( \underset{\_}{A} \right)}{\lambda_{\min}\left( \underset{\_}{A} \right)}};$however this property does not give any further insight intoperformance, since the matrix product, R _(x) ⊙R _(w), is biasedcompared to the ensemble covariance R _(x).

-   -   4. The analysis is for expected value, and assumes that the        {circumflex over (P)}_(y)(ψ) spectrum is the ideal P_(x)(ψ)        convolved with the function C_(w)(ψ). In practice, classical        power spectral density (PSD) techniques will average over the M        snapshots to reduce the variability in the estimate of the        spectrum.

Prototype Power Pattern

To assist in the normalized SINR loss analysis, a “prototype” windowfunction may be defined. A power pattern, C_(w)(ψ) may be defined by anideal bandlimited portion, which may be useful for identifying theimpacts of mainlobe width and a constant offset portion, which may beuseful for identifying the impacts of sidelobes or spectral leakage.Given definition for these regions, the analysis is more straightforwardand the two factors may be considered individually. Classically definedwindows incorporate both features together, in some trade-off related totheir design, making the analysis of individual window functions lessinsightful. The prototype window is defined herein using subscripts asfollows: lb denotes local bias and relates to the main lobe region, bbdenotes broadband bias and corresponds to the sidelobe levels. C_(w)(ψ)is a periodic function with a period 2π, and is defined explicitly overthe region of support, |ψ|≦π, asC _(w)(ψ)=C _(lb)(ψ)+C _(bb)(ψ)  (3.75)

The mainlobe is defined by an ideal bandlimited function

$\begin{matrix}{{C_{l\; b}(\psi)} = \left\{ \begin{matrix}A_{l\; b} & {{\psi } \leq \psi_{l\; b}} \\0 & {{\psi } > \psi_{l\; b}}\end{matrix} \right.} & (3.76)\end{matrix}$

and the sidelobe level is defined by a constantC _(bb)(ψ)=A _(bb)|ψ|≦π  (3.77)

The scale factors, A_(lb) and A_(bb), may be constrained according tothe normalization

$\begin{matrix}{{\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{C_{w}(\psi)}{\mathbb{d}\psi}}}} = 1} & (3.78)\end{matrix}$

By specifying two of the three parameters, {A_(lb), A_(bb), ψ_(lb)},e.g., the latter two, and varying them separately, their respectiveinfluence on performance can be examined. This provides a general feelfor behavior of normalized SINR loss using classical PSD techniques. Asexplained above, E{{circumflex over (R)} _(y)}=R _(x)

R _(w). Thus, it is sufficient to find R _(w) for the prototype window.For the definition in Eqn. (3.75), there is a closed form solution basedon the parameters. Starting with the Fourier transform relationship

$\begin{matrix}{{\rho_{w}\lbrack n\rbrack} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{C_{w}(\psi)}{\mathbb{e}}^{{j\psi}\; n}{\mathbb{d}\psi}}}}} & (3.79)\end{matrix}$

substituting in the definition for C_(w)(ψ) showing the explicit localand broadband bias terms

$\begin{matrix}{{\rho_{w}\lbrack n\rbrack} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{\left\lbrack {{C_{l\; b}(\psi)} + {C_{bb}(\psi)}} \right\rbrack{\mathbb{e}}^{j\;\psi\; n}{\mathbb{d}\psi}}}}} & (3.80)\end{matrix}$

The solution includes two terms, ρ_(w)[n]=ρ_(lb)[n]+ρ_(bb)[n]. The firstterm, ρ_(lb)[n], is:

$\begin{matrix}{{\rho_{l\; b}\lbrack n\rbrack} = {{\frac{1}{2\;\pi}{\int_{- \pi}^{\pi}{{C_{l\; b}(\psi)}{\mathbb{e}}^{{j\psi}\; n}{\mathbb{d}\psi}}}} = {{A_{l\; b}\left( {\psi_{l\; b}/\pi} \right)}\sin\;{c\left( {\left\lbrack {\psi_{l\; b}/\pi} \right\rbrack n} \right)}}}} & (3.81)\end{matrix}$

where sin c (x) ≡sin (πx)/(πx). The second term, ρ_(bb)[n], is

$\begin{matrix}{{\rho_{bb}\lbrack n\rbrack} = {{\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{C_{bb}(\psi)}{\mathbb{e}}^{{j\psi}\; n}{\mathbb{d}\psi}}}} = {A_{bb}{\delta\lbrack n\rbrack}}}} & (3.82)\end{matrix}$

The complete autocorrelation sequence, ρ_(w)[n], is the sum of the twoρ_(w) [n]=A _(lb)(ψ_(lb)/π)sin c([ψ_(lb) /π]m)+A _(bb) δ[m]  (3.83)

Expressed in matrix form, where R _(lb,bb)=((ρ_(lb,bb)[r−c]))_(r,c):R _(w) =R _(lb) +R _(bb)  (3.84)

The resultant expected value of the estimated covariance matrix is thenE{{circumflex over (R)} _(y) }=R _(x)

[R _(lb) +R _(bb)]  (3.85)

Observe in Eqn. (3.85) that the broadband bias term, R _(bb), is adiagonal matrix with constant diagonal A_(bb). The overall affect of thebroadband bias term is the same as a diagonal loading in that itincreases the main diagonal of the estimated covariance matrix. This isaccomplished as a multiplicative, not additive, effect.

The prototype window described by Eqn. (3.75) is not realizable with anyarray of finite length. So it is not possible to determine a set ofcoefficients, w, that would result in a power pattern of this type. Thisdoes not prevent analysis, however, since it is assumed that theestimated spectra is available, {circumflex over(P)}_(y)(ψ)=C_(w)(ψ)≠P_(x,ψ)(ψ), and not how it was computed. Theapproach is to understand performance issues when using classicaltechniques, and to show where they work well and where they do not.These results can then guide development to extend the range ofsituations allowing useful operation.

Prototype Window Normalized SINR Loss

Eqn (3.84) may be used in the expression for SINR loss of Eqn. (3.74).The Hadamard product structure within the expression of Eqn. (3.74)prevents a more compact form. Strictly speaking, the value of ξ foundusing Eqn. (3.74) is not a random variable but a constant, since theexpected value of the estimated covariance, E{{circumflex over(R)}_(y)}, is being used. To understand the statistical behavior of ξ,one may insert the computational form of the estimated covariance givenin Eqn. (3.65) into Eqn. (3.74). This is also not easily reducible to amore compact form.

The conclusion is that for this type of structured covariance matrixestimation algorithm, the normalized SINR loss cannot be simplified intoan expression that does not involve the particular covariance, R _(x),for the problem. This is unlike sample covariance matrix methods, wherethe performance can be derived in closed form as a random variable andshown to be function of number of elements and snapshots only.Performance of this structured covariance method depends on the problem.Eqn. (3.74) can still be used to understand SINR loss performance, butthe scenarios analyzed may be specified.

By specifying the parameters of the prototypical window function([A_(lb) or A_(bb)], ψ_(lb)), for particular types of interferenceproblems, Eqn (3.74) may be used to predict performance. For example,the single plane wave interferer in uncorrelated noise case may beconsidered with a signal of interest not present. This case highlightssome of the features of the approach.

Single Plane Wave Interferer

The impact of using CSS for the single plane wave in uncorrelated noisecase may be analyzed using Eqn. (3.74). Normalized SINR loss may becomputed as a function of distance between the desired signal directionof arrival and the interferer location in ψ-space, Δψ=(ψ_(s)−ψ_(o)). Foreach Δψ the difference between the per element INR and a fixed sidelobelevel, A_(bb), may be varied. This may be done, for example, for threevalues of ψ_(lb), set to multiples of the mainlobe half width of thearray,

$\left( {1,2,3} \right){\frac{2\pi}{N}.}$By increasing ψ_(lb) an indication of performance when using widermainlobe windows or when operating below the design frequency may bedetermined.

Simulation results indicate that for CSS covariance matrix estimationbased on classical power spectral based methods, near optimal SINR lossperformance can be achieved when the interference has sufficientseparation from the desired signal, or in the event the interferer isnear the desired signal, its power per sensor element is substantiallybelow the sidelobe level of the window used.

Analysis of a simple single plane wave interferer in white noise showedthat the adaptive processor is relatively insensitive to estimationerror of INR. Performance is more effected by wavenumber estimationaccuracy in direct relation to the accurate placement of nulls, butacceptable performance is achievable in few or even one snapshot. The CRbound is proportional to M⁻¹ and N⁻², so a larger number of sensors isbeneficial. This is in contrast to the closed form performance of samplecovariance matrix techniques, where performance does not depend on thenumber of sensors, just the number interferers and snapshots (with useof diagonal loading). Covariance matrix estimates developed usingclassical power spectral estimation techniques were seen to be biased.The normalized SINR loss performance indicated that to broaden theconditions under which useful performance can be achieved an additionalstep of harmonic analysis may be used, including the detection andsubtraction of line components from the data.

Structured Covariance Estimation with Multi-Taper Spectral Estimation(MTSE)

The performance of adaptive beamforming using covariance matrixestimates based on the wavenumber spectrum using classical spectralestimation techniques, e.g., power spectral distribution (PSD), isdescribed above. The expected value of the estimated covariance whenusing PSD is biased, E {R _(PSD)}=R _(x)⊙R _(w), so in general adaptivebeamformers based on the CSS covariance do not converge to the optimalsolution. However, analysis of the SINR loss performance for the uniformlinear array case, as function of interferer to desired signal spacing,Δψ, window characteristics, and INR showed that performance is within afew tenths of a dB from optimal under some conditions.

To maintain good normalized SINR loss performance, the interferer shouldhave sufficient separation from the desired signal, in proportion to thewindow mainlobe width, Δψ≧2ψ_(lb)−3ψ_(lb), with an interferer to noiseratio such that (INR_(PE)+A_(bb))≦10 log 10(N) (dB). To continue toachieve good performance with smaller separation, the condition(INR_(PE)+A_(bb))<<0 (dB) may be used. Neither of these conditions canbe guaranteed in practice. These concerns motivate use of multi-taperspectral estimation (MTSE), e.g., Thompson's MTSE, in forming theestimate of the frequency-wavenumber spectra instead of the classicaltechniques. This is for at least three reasons:

-   -   1. First, while the equivalent MTSE mainlobe width is larger        than traditional window functions it provides two benefits: 1)        averaging increased number of uncorrelated estimates of power        within a particular band in wavenumber-space results in lower        variance locally 2) estimates of the eigencoefficients within a        particular band in wavenumber-space allow for direct formulation        of a detection problem for line components based on the        underlying model for the stationary process and the Cramér        spectral representation theorem. This eliminates the need to        design new detection procedures based on estimated spectrum        only, such as local noise floor averaging, etc., and can be done        prior to forming the covariance matrix.    -   2. The parameters, amplitude, frequency and phase of the        detected line components may be used to coherently cancel        (subtract) them from the data and the estimation process may be        repeated on the residual. This process is known as harmonic        analysis and can be carried out iteratively.    -   3. Combination of the various eigenspectra may be accomplished        via an adaptive weighting mechanism designed to minimize local        and broadband bias, which may result in an equivalent window        with very low sidelobe level, e.g., below −60 dB.

Covariance from Spatial Spectra (CSS) with MTSE

The discussion below provides a procedural outline for using MTSE withharmonic analysis to estimate the frequency-wavenumber spectra used toform an estimate of the covariance matrix at the output of the array.The discussion considers an N element uniform linear array, andextension to more general geometries is considered further below. Theprocess has a mixed spectrum, with K point source signals haveindependent, random complex-valued amplitudes,

${a_{k}(m)}\overset{d}{->}{{{CN}\left( {0,\sigma_{k}^{2}} \right)}.}$M snapsnots are available for processing. The snapshot model is

$\begin{matrix}{{{\underset{\_}{x}}_{m} = {{\sum\limits_{k = 1}^{K}{{\underset{\_}{v}}_{k}{a_{k}(m)}}} + {\underset{\_}{n}}_{b,m} + {\underset{\_}{n}}_{w,m}}},{{\underset{\_}{x}}_{m}\overset{d}{\rightarrow}{CN}_{N}}} & (4.1)\end{matrix}$

An estimate for the frequency-wavenumber spectra, for both process andsensor noise, {circumflex over (P)}_(x)(ψ), may be formed using MTSE andused to compute an estimate of the covariance.

$\begin{matrix}{{\hat{\underset{\_}{R}}}_{x} = {\left( {2\pi} \right)^{- 1}{\int_{- \pi}^{\pi}{{{\hat{P}}_{x}(\psi)}{{\underset{\_}{v}}_{\psi}(\psi)}{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\mathbb{d}\psi}}}}} & (4.2)\end{matrix}$

The integral in Eqn. (4.2) may be implemented using a numericalsummation

$\begin{matrix}{{\hat{\underset{\_}{R}}}_{x} = {\frac{\Delta\; k}{2\pi}{\sum\limits_{b = {{- L}/2}}^{{L/2} - 1}{{{\hat{P}}_{x}\left( {l\;\Delta\; k} \right)}{\underset{\_}{v}\left( {l\;\Delta\; k} \right)}{{\underset{\_}{v}}^{H}\left( {l\;\Delta\; k} \right)}}}}} & (4.3)\end{matrix}$

Alternatively, for a uniform linear array one can take advantage of thefrequency-wavenumber spectrum and covariance being single dimension, andrelate them by the 1-D inverse Fourier transform

$\begin{matrix}{{{\hat{\underset{\_}{R}}}_{x}\lbrack n\rbrack} = {\frac{1}{2\pi}{\sum\limits_{l = 0}^{L - 1}{{{\hat{P}}_{x}\left( {l\;\Delta\; k} \right)}{\mathbb{e}}^{j\; l\;\Delta\;{kn}}}}}} & (4.4)\end{matrix}$

and populate the covariance matrix as {circumflex over (R)}_(x)=(({circumflex over (R)}_(x)[r−c]))_(r,c). Fast Fourier transformtechniques may be a particularly efficient implementation when the valueof L permits their use.

Number of Tapers, D

To exploit the efficiencies of FFT based processing when using MTSE toestimate the spectrum, {circumflex over (P)}_(x)(ψ), the spectralestimation may maintain a fixed resolution in normalized wavenumberspace, ψ=−k_(z)d. This allows a single set of tapers to be used. A widthof the analysis region, W, may be chosen. Typically choices are NW=1.5,2.0, 2.5. The number of significant tapers supported by a particularwidth, W, is D=2 NW−1. The case of NW=1 results in a single usable taperand reverts to a standard classical PSD method. No harmonic analysis maybe used in this scenario. As a practical matter, D=2 NW may be selected,resulting in a lower concentration for the last taper. This may increasethe number of basis vectors used in the harmonic analysis detectionstatistic. The D tapers may be designed according to an appropriateeigenvalue or generalized eigenvalue problem. For the ULA case, theresultant tapers are a discrete prolate spheroidal sequences.q _(d)(n)=dpss(N,NW),d=1, . . . ,D  (4.5)

FFT and Zero-Padding

Snapshot data may be windowed and fast Fourier transformed to producethe MTSE eigencoefficients

$\begin{matrix}{{y_{m}^{(d)}(l)} = {\sum\limits_{n = 0}^{N - 1}{{x_{m}(n)}{q_{d}(n)}{\exp\left( {{- {j2}}\;\pi\;{\ln/N_{FFT}}} \right)}}}} & (4.6)\end{matrix}$

for the set of points l=0, . . . ,N_(FFT)−1. The FFT size, N_(FFT), maybe specified in Eqn. (4.6) independently from the number of arrayelements, N. The nominal set of points would be N_(FFT)=N, with agreater number of points, N_(FFT)>N, generated using the zero-paddingtechnique (N_(FFT)<N is possible using polyphase techniques but notlikely a case of interest here). The zero-pad operation may be usefulfor several reasons. First, the detection process within harmonicanalysis performs a subtraction of discrete harmonic (point source)components in the data. This is done by estimating the unknown linecomponent parameters: wavenumber, amplitude, and phase. Many algorithmsexist for estimating parameters of sinusoids in noise. Assuming thatmultiple interferers are sufficiently separated, these parameters may beconveniently estimated optimally using FFT techniques. The precision towhich this can be accomplished, without additional techniques, such ascurve fitting between FFT bins, is directly proportional to the finenessof the FFT spacing in wavenumber space. The zero-padding operation is anefficient method for increasing this fineness. The zero-padding is alsouseful in the smoothing operation referred to as free parameterexpansion, FPE.

Discrete Line Component Processing (Harmonic Analysis)

Harmonic analysis algorithm operates on the eigencoefficients,Y_(m,d)(l), to determine the presence of discrete line components. Adetection statistic is computed as a ratio of the power in the linecomponent subspace to the power outside that subspace in the region[f_(o)−W<f≦f_(o)+W].

$\begin{matrix}{{F\left( f_{o} \right)} = {\frac{\sum\limits_{m = 1}^{M}{{{\underset{\_}{y}}_{m}^{H}\left( f_{o} \right)}{\underset{\_}{P}}_{q}{{\underset{\_}{y}}_{m}\left( f_{o} \right)}}}{\sum\limits_{m = 1}^{M}{{{\underset{\_}{y}}_{m}^{H}\left( f_{o} \right)}{\underset{\_}{P}}_{q}^{\bot}{{\underset{\_}{y}}_{m}\left( f_{o} \right)}}}\overset{H_{1}}{\underset{H_{o}}{\gtrless}}{\Upsilon\;{TH}}}} & (4.7)\end{matrix}$

The choice of threshold, γTH, can be determined using a Neyman-Pearsoncriteria assuming Gaussian noise statistics. Practically, it may also beuseful to define a minimum limit allowable for detection, e.g., γ_(min)^((dB))=10 log₁₀(γ_(min)), such thatγTH=max(γ_(NP),γ_(min))  (4.8))

with γ_(min) ^((dB))=3 dB typically. This reduces excessive falsedetections and has minimal impact on overall performance, since theharmonic analysis is used primarily to eliminate high powered, not lowpowered, discrete interference. This test is valid for a single linecomponent present within the analysis region, [f_(o)−W<f≦f_(o)+W]. Ifthe interference environment is dense with respect to the arrayresolution, additional tests such as a double F line test may beappropriate.

For a detected line component, it may be assumed that the wavenumberremains constant across all snapshots. For the k^(th) line component thewavenumber, ψ_(k), is used to estimate the remaining parameters persnapshot and form the overall covariance matrix. With sufficientzero-padding, ψ_(k), can be estimated as

$\begin{matrix}{{\hat{\psi}}_{k} = {{\underset{\psi}{argmax}{F(\psi)}} > {\Upsilon\;{TH}}}} & (4.9)\end{matrix}$

The remaining parameters may be estimated per snapshot using matchedfilter techniques. Defining the reference sinusoid waveforms _(ref,k)=((e ^(jψ) ^(k) ^(n)))_(n)  (4.10)

the complex-amplitude and interferer power are estimated as

$\begin{matrix}{{{{\hat{a}}_{k}(m)} = {\frac{1}{N}{\underset{\_}{s}}_{{ref},k}^{H}{\underset{\_}{x}}_{m}}},{{\hat{\sigma}}_{k}^{2} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{{\hat{a}}_{k}(m)}}^{2}}}}} & (4.11)\end{matrix}$

While not required for the processing, the interferer to noise ratio foreach detected discrete component may be estimated. This is useful whengenerating a composite spectrum for visualization. A composite spectrummay be generated based on the estimated continuous background spectrumof the residual snapshot data (post harmonic analysis), with numericalinsertion of the discrete components. The insertion technique may useknowledge of the INR to properly represent the uncertainty of aparticular estimate.

$\begin{matrix}{{I\hat{N}R_{k}} = \frac{\sum\limits_{m = 1}^{M}{{{\hat{a}}_{k}(m)}}^{2}}{\left( {D - 1} \right)^{- 1}{\sum\limits_{m = 1}^{M}{{{\underset{\_}{y}}_{m}^{H}\left( {\hat{\psi}}_{k} \right)}{\underset{\_}{P}}_{q}^{\bot}{{\underset{\_}{y}}_{m}\left( {\hat{\psi}}_{k} \right)}}}}} & (4.12)\end{matrix}$

Influence of the K detected components may be removed from the snapshotdata to produce the residual data snapshots, x _(b,m). This can beaccomplished with one of two methods.

Method 1. Coherent Subtraction

$\begin{matrix}{{\underset{\_}{x}}_{b,m} = {{\underset{\_}{x}}_{m} - {\sum\limits_{k = 1}^{K}{{{\hat{a}}_{k}(m)}{\underset{\_}{s}}_{{ref},k}}}}} & (4.13)\end{matrix}$

Method 2. Null Projection

$\begin{matrix}{{\underset{\_}{P}}_{\underset{k}{\bot}} = {\underset{\_}{I} - {{{\underset{\_}{s}}_{{ref},k}\left( {{\underset{\_}{s}}_{{ref},k}^{H}{\underset{\_}{s}}_{{ref},k}} \right)}^{- 1}{\underset{\_}{s}}_{{ref},k}^{H}}}} & (4.14) \\{{\underset{\_}{P}}_{k}^{\bot} = {\prod\limits_{k = 1}^{K}{\underset{\_}{P}}_{k}^{\bot}}} & (4.15)\end{matrix}$x _(b,m) =P _(K) ^(⊥) x _(m)  (4.16)

Continuous Background Spectrum

Once harmonic analysis is complete the residual snapshot data may beused to compute a final smooth, continuous background spectrum. Thesnapshot data may be windowed and fast Fourier transformed to produceeigencoefficients

$\begin{matrix}{{y_{b,m}^{(d)}(l)} = {\sum\limits_{n = 0}^{N - 1}{{x_{b,m}(n)}{q_{d}(n)}{\exp\left( {{- {j2\pi}}\;{\ln/N_{FFT}}} \right)}}}} & (4.17)\end{matrix}$

The eigencoefficients may be used to produce the individualeigenspectra.

$\begin{matrix}{{{\hat{P}}_{b}^{(d)}(l)} = {{y_{b,m}^{(d)}(l)}}^{2}} & (4.18)\end{matrix}$

The individual eigenspectra may be linearly combined according to a setof weights

$\begin{matrix}{{{\hat{P}}_{x,b}(l)} = {\sum\limits_{d = 1}^{D}{{h_{d}(l)}{{\hat{P}}_{b}^{(d)}(l)}}}} & (4.19)\end{matrix}$

The weights, h_(d)(k), may be fixed (e.g., for an underlying whitespectrum) or may be determined adaptively.

Covariance Matrix Computation

The final estimate of the covariance matrix may be formed using linecomponent and continuous background spectrum products.{circumflex over (R)} _(a)=diag(σ₁ ²,σ₂ ², . . . ,σ_(K) ²)  (4.20){circumflex over (V)}=[ s _(ref,1) ,s _(ref,2) , . . . ,s_(ref,K)]  (4.21)

$\begin{matrix}{{\hat{\underset{\_}{R}}}_{MTSE} = {{{\underset{\_}{\hat{V}\hat{R}}}_{a}{\underset{\_}{\hat{V}}}^{H}} + {\left( {2\pi} \right)^{- 1}{\int_{- \pi}^{\pi}{{{\hat{P}}_{x,b}(\psi)}{{\underset{\_}{v}}_{\psi}(\psi)}{{\underset{\_}{v}}_{\psi}^{H}(\psi)}{\mathbb{d}\psi}}}}}} & (4.22)\end{matrix}$

Composite Spectrum Generation

In a particular embodiment, the wavenumber spectrum, {circumflex over(P)}_(x) (ψ), may be visualized directly in addition to forming thecovariance matrix. The smooth continuous component is the direct outputfrom the MTSE processing of the residual snapshot data, x _(b,m),yielding {circumflex over (P)}_(x,b)(ψ). The discrete componentspreviously estimated and subtracted are then added back into thenumerical spectrum. Two points that may be considered in this contextare: 1) numerical power should conserved and 2) better resolution orperformance than is available with the data should not be implied.

Unbiased Spectral Estimate in White Noise

In this type of spectral estimation, each point in the spectral estimateis scaled such that it is an unbiased estimate of the noise power for awhite noise input. In terms of classical PSD techniques, this impliesthat the window function has been scaled such that w ^(H) w=1.0. MTSEtapers may be scaled in accordance with this approach. Plane wave ordiscrete sinusoidal components experience a processing gain due to thecoherent gain of the window function. The maximum gain may be obtainedusing

${\underset{\_}{w} = {\frac{1}{\sqrt{N}}\underset{\_}{1}}},$and equates to 10 log 10(N) dB. Thus, a snapshot,

${{\underset{\_}{x}}_{m}\overset{d}{->}{{CN}_{N}\left( {0,{\sigma^{2}\underset{\_}{I}}} \right)}},$will produce an expected value at any power spectral estimate, E{{circumflex over (P)} (ψ)}=σ2, while a snapshot corresponding to aplane wave component, x _(m)=A_(o) v _(ψ)(ψ_(o)) will produce a value ofE {{circumflex over (P)}(ψ_(o))}≦NA_(o) ². The coherent gain of the MTSEtapers can be found by computing

$\begin{matrix}{{CG}_{MTSE} = {\frac{1}{D}{\sum\limits_{d = 1}^{D}{{\sum\limits_{n = 0}^{N - 1}{q_{d}(n)}}}^{2}}}} & (4.23)\end{matrix}$

Because the choice of coherent gain is somewhat arbitrary for there-inserted spectral component, it may be convenient to use the maximumtheoretical processing gain, CG=10 log₁₀ (N), such that the compositespectrum looks similar, on a relative scale, to that obtained using MVDRtechniques.

$\begin{matrix}{{{\hat{P}}_{dB}\left( \psi_{k} \right)} = {{10{\log_{0}\left( {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{A_{k}(m)}}^{2}}} \right)}} + {10{\log_{10}(N)}}}} & (4.24)\end{matrix}$

Normalized SINR loss performance was assessed via simulation for anumber of interference and noise scenarios. Performance was seen toconverge with very few snapshots to near optimal, and in many cases withfewer snapshots than interferers. These simulation results aresignificantly better than what is achievable using diagonal loading orcomparable reduced rank techniques. The simulations considered linecomponent, spatially spread, and mixed spectra conditions. It wasobserved that performance could be improved by increasing the estimationaccuracy of the harmonic analysis step to provide better cancellation ofthe line components within the data.

Correlated Signal and Interference

One of the properties of a wide sense stationary narrowband space-timeprocess (as considered above) is that it can be represented as a sum ofuncorrelated plane waves, distributed across all directions of arrivalto the array. This model does not describe every situation that may beencountered. An example of another case of interest is when there iscorrelation between two or more plane wave components observed by thearray. This may occur in situations of multipath or smart jamming. Underthese conditions the covariance is function of absolute position, notrelative, and the process is not wide sense stationary. Failure toaccount for the correlation within the data can lead to signalcancellation, and an overall loss in output SNR.

In a particular embodiment, the snapshot data may include both thesignal of interest and interference correlated with it. Wavenumber orspatial spectra provide no correlation information. Covariance fromspatial spectrum (CSS) may provide a level of robustness against theeffects of correlated signal and interference on performance. CSS may bebiased in two ways. A first bias is in a manner similar to the biasdiscussed above for the wide sense stationary process case due to themethod of spectral estimation. The second bias is specific to thecorrelation within the data. Performance may be assessed in comparisonto CSS operating on uncorrelated data, as well with the effective SINRmetric, an appropriate measure of output SINR in the correlated signaland interference scenario. Using simulations, the bias attributable tothe correlation component was found to have negligible impact onperformance.

Covariance for Correlated Signals

In the case where the point source signals are correlated the space-timeprocess is not spatially stationary. Referring back to the Cramérspectral representation of the stationary space-time process, thecorrelation violates the condition that disjoint regions in wavenumberspace be uncorrelated. In terms of the covariance, the effect is seen asa dependence on absolute as well as relative position. This is visiblewhen examining the covariance matrix based on the model for the snapshotdata, x _(m).

$\begin{matrix}{{\underset{\_}{x}}_{m} = {{{\sum\limits_{k = 1}^{K}{{\underset{\_}{v}}_{k}{a_{k}(m)}}} + {\underset{\_}{n}}_{m}} = {{\underset{\_}{Va}}_{m} + {\underset{\_}{n}}_{m}}}} & (5.1)\end{matrix}$

where the background noise component and sensor noise component arecombined together,

${{\underset{\_}{n}}_{m} = {{\underset{\_}{n}}_{b,m} + {\underset{\_}{n}}_{w,m}}},{{{\underset{\_}{n}}_{m}\overset{d}{\longrightarrow}{{CN}_{N}\left( {\underset{\_}{0},{\underset{\_}{R}}_{n}} \right)}}.}$The covariance matrix isE{x _(m) x _(m) ^(H) }=R _(x) =VR _(a) V ^(H) +R _(n)  (5.2)where R _(a)=E {a _(m) a _(m) ^(H)}. At least in this context, R _(a) isnot a diagonal matrix. The off diagonal terms represent thecross-correlation between the plane waves. R _(a) can be expressed as acombination of a diagonal matrix and an off-diagonal matrix.R _(a) =R _(a,U) +R _(a,C)  (5.3)

The subscript U is used to reinforce that the diagonal matrix relates touncorrelated plane waves, while the subscript C is used to reinforcethat the off-diagonal matrix corresponds to the terms representing thecorrelation. Thus,R _(x) =VR _(a,U) V ^(H) +R _(n) +R _(a,C) V ^(H)  (5.4)

The first two terms of Eqn. (5.4) correspond to the covariance for thestationary process model. This portion of the overall covariance as isreferred to herein as:R _(x,U) =VR _(a,U) V ^(H) +R _(n)  (5.5)

The remaining portion, due to the off-diagonal entries of the matrix R_(a), is referred to in a similar fashion.R _(x,C) =VR _(a,C) V ^(H)  (5.6)

Expected Value

To proceed in analyzing the expected value of the covariance matrixestimate, CSS with classical power spectral estimation may be used. Asexplained above, for a stationary process:

$\begin{matrix}{{{\underset{\_}{x}}_{m} = {{\sum\limits_{k = 1}^{K}{{\underset{\_}{v}}_{k}{a_{k}(m)}}} + {\underset{\_}{n}}_{m}}},{{\underset{\_}{x}}_{m}\overset{d}{\longrightarrow}{{CN}_{N}\left( {\underset{\_}{0},{{{\underset{\_}{VR}}_{a,U}{\underset{\_}{V}}^{H}} + {\underset{\_}{R}}_{n}}} \right)}}} & (5.7)\end{matrix}$

where the subscript U is used to reinforce that the process isstationary and that R _(a,U) is a diagonal matrix. The CSS covariancematrix estimate may be found from the windowed (tapered) snapshots, y_(m)=x _(m)⊙w, as

$\begin{matrix}{{\underset{\_}{\hat{R}}}_{y} = {{{DSR}\left( {\underset{\_}{R}}_{w,{SCM}} \right)} = {{DSR}\left( {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{\underset{\_}{y}}_{m}{\underset{\_}{y}}_{m}^{H}}}} \right)}}} & (5.8)\end{matrix}$

with expected value E {{circumflex over (R)} _(y)}=R _(x,U)

R _(w). Alternatively, the windowed snapshot model for this case may bewritten as:

$\begin{matrix}\begin{matrix}{{\underset{\_}{y}}_{m} = {{\sum\limits_{k = 1}^{K}{\left( {{\underset{\_}{v}}_{k}\bullet\underset{\_}{w}} \right){a_{k}(m)}}} + \left( {{\underset{\_}{n}}_{m}\bullet\underset{\_}{w}} \right)}} \\{= {{\sum\limits_{k = 1}^{K}{{\overset{\_}{\underset{\_}{v}}}_{k}{a_{k}(m)}}} + {\overset{\_}{\underset{\_}{n}}}_{m}}} \\{= {{\underset{\_}{\overset{\_}{V}a}}_{m} + {\overset{\_}{\underset{\_}{n}}}_{m}}}\end{matrix} & (5.9)\end{matrix}$

so that the tapered snapshots are distributed as

$\begin{matrix}{{\underset{\_}{y}}_{m}\overset{d}{\longrightarrow}{{CN}_{N}\left( {\underset{\_}{0},{{{\underset{\_}{\overset{\_}{V}R}}_{a,U}{\overset{\_}{\underset{\_}{V}}}^{H}} + {\overset{\_}{\underset{\_}{R}}}_{n}}} \right)}} & (5.10)\end{matrix}$

From Eqn. (5.8), DSR linearity properties and Eqn. (5.10):

E {{circumflex over (R)} _(y)}=DSR E {y _(m) y _(m) ^(H)})=DSR( VR_(a,U) V ^(H)+ R _(n))

The two expressions for the expectation are equivalent. Thus,R _(x,U)

R _(w)=DSR( VR _(a,U) V ^(H)+ R _(n))  (5.11)

The plane wave components may be correlated, so that

$\begin{matrix}{{\underset{\_}{y}}_{m}\overset{d}{\longrightarrow}{{CN}_{N}\left( {\underset{\_}{0},{{{\overset{\_}{\underset{\_}{V}}\left\lbrack {{\underset{\_}{R}}_{a,U} + {\underset{\_}{R}}_{a,C}} \right\rbrack}{\overset{\_}{\underset{\_}{V}}}^{H}} + {\overset{\_}{\underset{\_}{R}}}_{n}}} \right)}} & (5.12)\end{matrix}$

Eqn. (5.8) may be used, so that the expected value of the estimatedcovariance isE{{circumflex over (R)} _(y)}=DSR( V [ R _(a,U) +R _(a,C) V ^(H)+ R_(n)  (5.13)

Using the DSR properties, the expected value includes two termsE{{circumflex over (R)} _(y)}=DSR( VR _(a,U) V _(H)+ R _(n))+DSR( VR_(a,C) V _(H))  (5.14)

From Eqn. (5.11), the first term is the CSS with classical spectralestimation covariance as if the process where in fact stationary, so thefinal result isE{{circumflex over (R)} _(y}) =R _(x,U)

R _(w)+DSR( VR _(a,C) V ^(H)  (5.15)

Eqn. (5.15) implies that for a correlated signal and interferenceproblem, the CSS with classical spectral estimation technique produces acovariance matrix estimate that is an estimate of the covariance as ifthe process were uncorrelated, R _(x,U)⊙R _(w), with an additional biasterm DSR ( VR _(a,C) V ^(H)). In terms of addressing signalcancellation, the first term has eliminated the correlation in theestimated covariance. Thus, an adaptive beamformer based on this part ofthe estimate would be expected to have a performance consistent with thediscussion above, even if the process itself has correlated components.The impact of the remaining bias term on performance is discussedfurther below.

CSS Performance with Correlated Signal and Interference

Correlated signal and interference introduces the potential for signalcancellation for some adaptive beamforming algorithms. A minimumvariance distortionless response (MVDR) beamformer may be derived to beoptimal when observing uncorrelated noise and interference only, aspatially stationary space-time process. The MVDR approach can beextended when both the desired signal and interference are present inthe snapshot data, the so called minimum power distortionless response(MPDR) beamformer. MPDR attempts to minimize output power whileconstrained to be distortionless in the direction of the desired signal,w _(MPDR) ∝R ¹ s. This allows the desired signal through due to thedistortionless constraint, but in the event of correlated interferencethe processor may use the interferer to destructively cancel the desiredsignal in the overall attempt to minimize output power.

Relative Contribution of Correlated and Uncorrelated Components

Questions that may be of interest for adaptive beamformers based uponthe covariance from spatial spectrum estimate, referred to as{circumflex over (R)} _(CSS), include: 1) does {circumflex over (R)}_(CSS) convey any information regarding the correlation component in thedata if it exists, and 2) how well do beamformers based upon {circumflexover (R)}_(CSS) perform compared to an adaptive processor using acovariance for the data where there is no correlation present? Asexplained above, the expected value of {circumflex over (R)} _(CSS)contains the covariance if the data were uncorrelated plus an additionalterm.

$\begin{matrix}\begin{matrix}{{E\left\{ {\hat{\underset{\_}{R}}}_{CSS} \right\}} = {{{\underset{\_}{R}}_{x,U}\mspace{11mu}\bullet\mspace{11mu}{\underset{\_}{R}}_{w}} + {{DSR}\left( {\underset{\_}{\overset{\_}{V}}\;{\underset{\_}{R}}_{a,C}{\underset{\_}{\overset{\_}{V}}}^{H}} \right)}}} \\{= {{\underset{\_}{R}}_{{CSS},U} + {\underset{\_}{R}}_{{CSS},C}}}\end{matrix} & (5.16)\end{matrix}$

The influence of the second temi can be considered a bias. Byconsidering the case of two interferers with no noise component,explicit expressions for the bias term show that for redundancyaveraging the bias is not guaranteed to go to zero even if the arraylength is extended infinitely. Using the ratio of the Froebenius normsquared, ∥∥_(F) ², of the bias component to the unbiased component, asarray length is extended to infinity, the ratio goes to zero. Thisimplies that while the bias term exists, its impact as measured by therelative power indicated by ∥∥_(F) ² may become insignificant as thearray length increases.

A similar line of analysis may be followed using a two tone scenariowith no noise. The signals may have angles of arrival denoted as ψ₁,ψ₂with corresponding array manifold response vectors, v ₁, v ₂, andrespective variances σ₁ ², σ₂ ². The correlation between the two may bedescribed by a magnitude and phase as ζ=E {a₁ (m)a₂*(m)}=A_(ζ)e^(j)

^(ζ). The ensemble covariance may include correlated and uncorrelatedterms.

$\begin{matrix}\begin{matrix}{\underset{\_}{R} = {{\underset{\_}{V}{\underset{\_}{R}}_{a,U}{\underset{\_}{V}}^{H}} + {\underset{\_}{V}{\underset{\_}{R}}_{a,C}{\underset{\_}{V}}^{H}}}} \\{= {{\underset{\_}{R}}_{U} + {\underset{\_}{R}}_{C}}}\end{matrix} & (5.17)\end{matrix}$

Of interest is where the relative ratio of ∥∥_(F) ² shows that the biasterm is insignificant within E {{circumflex over (R)} _(CSS)}, comparedto the same ratio for the ensemble covariance,

$\begin{matrix}{\frac{{{\hat{\underset{\_}{R}}}_{{CSS},C}}_{F}^{2}}{{{\hat{\underset{\_}{R}}}_{{CSS},U}}_{F}^{2}}\bullet\frac{{{\underset{\_}{R}}_{C}}_{F}^{2}}{{{\underset{\_}{R}}_{U}}_{F}^{2}}} & (5.18)\end{matrix}$

When Eqn. (5.18) is valid, it indicates that CSS has substantiallydiminished the contributed of the correlated component, as measuredusing ∥∥_(F) ². Consider the problem for an N element uniform lineararray. The uncorrelated and correlated components for the ensemblecovariance are

$\begin{matrix}{\mspace{79mu}{{{\underset{\_}{R}}_{U}}_{F}^{2} = {{N^{2}\left( {\sigma_{1}^{4} + \sigma_{2}^{4}} \right)} + {2\sigma_{1}^{2}\sigma_{2}^{2}{{{\underset{\_}{v}}_{1}^{H}{\underset{\_}{v}}_{2}}}^{2}}}}} & (5.19) \\{{{\underset{\_}{R}}_{C}}_{F}^{2} = {2\sigma_{1}^{2}\sigma_{2}^{2}{A_{\zeta}^{2}\left\lbrack {N^{2} + {{\cos\left( {{2{\bullet\zeta}} + {{\Delta\psi}\left\lbrack {N - 1} \right\rbrack}} \right)}\frac{\sin^{2}\left( {{\Delta\psi}\;{N/2}} \right)}{\sin^{2}\left( {{\Delta\psi}/2} \right)}}} \right\rbrack}}} & (5.20)\end{matrix}$

where Δψ=ψ₂−ψ₁. The expressions for the components of E {{circumflexover (R)} _(CSS)} use the following property. For a Hermitian, Toeplitzmatrix,

$\begin{matrix}{\underset{\_}{A} = \begin{pmatrix}{\rho\lbrack 0\rbrack} & {\rho^{*}\lbrack 0\rbrack} & {\rho^{*}\lbrack 2\rbrack} & \ldots \\{\rho\lbrack 1\rbrack} & {\rho\lbrack 0\rbrack} & {\rho^{*}\lbrack 1\rbrack} & \; \\{\rho\lbrack 2\rbrack} & {\rho\lbrack 1\rbrack} & {\rho\lbrack 0\rbrack} & \ddots \\\vdots & \; & \ddots & \ddots\end{pmatrix}} & (5.21)\end{matrix}$

The Froebenius norm squared of A is

$\begin{matrix}{{\underset{\_}{A}}_{F}^{2} = {\sum\limits_{n = {- {({N - 1})}}}^{N - 1}{\left( {N - {n}} \right){{\rho\lbrack n\rbrack}}^{2}}}} & (5.22)\end{matrix}$

This yields a simplified expression for the ∥∥_(F) ² of the uncorrelatedcomponent of E {{circumflex over (R)} _(CSS)}.

$\begin{matrix}{{{\hat{\underset{\_}{R}}}_{{CSS},U}}_{F}^{2} = {\sum\limits_{n = {- {({N - 1})}}}^{({N - 1})}{\left( {N - {n}} \right){{{\rho_{w}\lbrack n\rbrack}}^{2}\left\lbrack {\left( {\sigma_{1}^{4} + \sigma_{2}^{4}} \right) + {2\sigma_{1}^{2}\sigma_{2}^{2}{\cos\left( {{\Delta\psi}\; n} \right)}}} \right\rbrack}}}} & (5.23)\end{matrix}$

where ρ_(w)[n] is the sample autocorrelation of the taper used, w. Thecorrelated component cannot be similarly reduced because of thesummation term in the sample autocorrelation, ρ_(CSS,C) [n], with asimplest expression for arbitrary w given as

$\begin{matrix}{{{\hat{\underset{\_}{R}}}_{{CSS},C}}_{F}^{2} = {\sum\limits_{n = {- {({N - 1})}}}^{N - 1}{\left( {N - {n}} \right){{\rho_{{CSS},C}\lbrack n\rbrack}}^{2}}}} & (5.24)\end{matrix}$

where

$\begin{matrix}{{{\rho_{{CSS},C}\lbrack n\rbrack} = {\sum\limits_{r = n}^{N - 1}\;{\sigma_{1}\sigma_{2}{w(r)}{{w^{*}\left( {r - n} \right)}\left\lbrack {{\zeta\mathbb{e}}^{{{j\psi}_{2}r} - {\psi_{1}{\lbrack{r - n}\rbrack}}} + {\zeta^{*}{\mathbb{e}}^{{{j\psi}_{1}r} - {\psi_{2}{\lbrack{r - n}\rbrack}}}}} \right\rbrack}}}}\mspace{79mu}{{n \geq 0},\;{{\rho_{{CSS},C}\left\lbrack {- n} \right\rbrack} = {\rho_{{CSS},C}^{*}\lbrack n\rbrack}}}} & (5.25) \\{{{\rho_{{CSS},C}\lbrack n\rbrack}}^{2} = {A_{\zeta}^{2}\sigma_{1}^{2}\sigma_{2}^{2}{{\sum\limits_{r = n}^{N - 1}\;{{w(r)}{{w^{*}\left( {r - n} \right)}\left\lbrack {{\mathbb{e}}^{{{j\psi}_{2}r} - {\psi_{1}{\lbrack{r - n}\rbrack}} + {\bullet\zeta}} + {\rho^{*}{\mathbb{e}}^{{{j\psi}_{1}r} - {\psi_{2}{\lbrack{r - n}\rbrack}} - {\bullet\zeta}}}} \right\rbrack}}}}^{2}}} & (5.26)\end{matrix}$

Normalized SINR Loss w.r.t. Uncorrelated Ensemble Covariance

Normalized SINR loss provides a useful performance metric for practicalarray lengths. Because MVDR is not designed for the correlated signaland interference case, normalized SINR loss calculations based on theensemble covariance containing the correlation may be inappropriate. Forexample, they do not predict the detrimental effects of the signalcancellation described earlier. The normalized SINR loss may be used tounderstand performance, and the ensemble covariance may be used for theuncorrelated data scenario, R _(x,U), as the reference covariance. Thiseffectively compares the performance of CSS with the optimal beamformeras if there were no correlation contained in the data. The subscript R_(C) is used to indicate that this is the normalized SINR loss forcorrelated data case.

$\begin{matrix}{\xi_{{\underset{\_}{R}}_{c}} = {\frac{{{\underset{\_}{s}}^{H}\left( {{\hat{\underset{\_}{R}}}_{{CSS},U} + {\hat{\underset{\_}{R}}}_{{CSS},C}} \right)}^{- 1}\underset{\_}{s}}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{x,U}^{- 1}\underset{\_}{s}} \times \frac{{{\underset{\_}{s}}^{H}\left( {{\hat{\underset{\_}{R}}}_{{CSS},U} + {\hat{\underset{\_}{R}}}_{{CSS},C}} \right)}^{- 1}\underset{\_}{s}}{{{\underset{\_}{s}}^{H}\left( {{\hat{\underset{\_}{R}}}_{{CSS},U} + {\hat{\underset{\_}{R}}}_{{CSS},C}} \right)}^{- 1}{{\underset{\_}{R}}_{x,U}\left( {{\hat{\underset{\_}{R}}}_{{CSS},U} + {\hat{\underset{\_}{R}}}_{{CSS},C}} \right)}^{- 1}\underset{\_}{s}}}} & (5.28)\end{matrix}$

As a second measure of performance in addition to ξ _(R) _(C) , which isan absolute measure, the relative performance compared to CSS when thereis no correlation within the data may be of interest. This performanceis referred to as ξ _(R) _(C) , where

$\begin{matrix}{\xi_{{\underset{\_}{R}}_{U}} = {\frac{{{\underset{\_}{s}}^{H}\left( {\hat{\underset{\_}{R}}}_{{CSS},U} \right)}^{- 1}\underset{\_}{s}}{{\underset{\_}{s}}^{H}{\underset{\_}{R}}_{x,U}^{- 1}\underset{\_}{s}} \cdot \frac{{{\underset{\_}{s}}^{H}\left( {\hat{\underset{\_}{R}}}_{{CSS},U} \right)}^{- 1}\underset{\_}{s}}{{{\underset{\_}{s}}^{H}\left( {\hat{\underset{\_}{R}}}_{{CSS},U} \right)}^{- 1}{{\underset{\_}{R}}_{x,U}\left( {\hat{\underset{\_}{R}}}_{{CSS},U} \right)}^{- 1}\underset{\_}{s}}}} & (5.29)\end{matrix}$

The difference between the performance measures of Eqn. (5.28) and Eqn.(5.29) isΔξ_(dB) =ξR _(C) _(,dB) −ξ _(R) _(U) _(,dB)   (5.30)

where ξ_(dB)=−10 log₁₀ξ.

Simulation was performed for a case of an N=32 uniform linear observingcorrelated signal and interference. The simulation determined apredicted normalized SINR loss (dB) of CSS with classical spectralestimation as a function of separation between the sources and thecorrelation angle with correlation magnitude fixed at ζ=1. Performancewas determined with respect to an optimal processor using an ensemblecovariance where the signal and interference were uncorrelated. Such aprocessor may be the best possible when a goal is to reduce the impactof the correlation in the data. The simulated performance was within 0.2dB almost everywhere, except when the sources were very close. Thus, thesimulations indicate that CSS may perform nearly identically whether thesignal and interference are correlated or not.

Results of the simulations indicate that while E{{circumflex over(R)}_(CSS)} has a bias component when correlated signal and interferenceare present, the impact of the bias component is negligible asperformance is nearly identical to the uncorrelated signal andinterference case. Covariance estimates based on frequency wavenumberspectrum have a decorrelating effect when signal and interference arecorrelated, and reduce the potential for signal cancellation tonegatively impact performance. These results indicate that CSStechniques should perform consistently in the presence of correlatedsignal and interference, regardless of correlation coefficient or angle.Also, the performance should be in line with CSS performance as if thedata were uncorrelated.

Redundancy Averaging

Redundancy averaging may be used to address a correlated signal andinterference problem. Redundancy averaging takes advantage of themultiple available estimates of the space-time correlation at a givenspatial lag by averaging them, and then generates a covariance matrixusing the averaged values. For the uniform line array, this amounts toreplacing diagonals in the sample covariance matrix with the averagediagonal values. For redundancy averaging, an averaged diagonal valuemay be defined for the n^(th) diagonal for the m^(th) snapshot, x_(m)=((x_(m)[n]))_(n)

$\begin{matrix}{{\rho_{{RA},m}\lbrack n\rbrack} = {\frac{1}{N - {n}}{\sum\limits_{\beta = 0}^{N - 1}\;{{x_{m}\lbrack\beta\rbrack}{x_{m}^{*}\left\lbrack {\beta - n} \right\rbrack}}}}} & (5.31)\end{matrix}$

where it is assumed that the sequence x_(m)[β]=0 for β<0,β≧N. The sampleautocorrelation for the each snapshot, x _(m), is

$\begin{matrix}{{\rho_{x,m}\lbrack n\rbrack} = {\sum\limits_{\beta = 0}^{N - 1}\;{{x_{m}\lbrack\beta\rbrack}{x_{m}^{*}\left\lbrack {\beta - n} \right\rbrack}}}} & (5.32)\end{matrix}$

Eqn. (5.31) applies a sample autocorrelation, ρ _(w) [n], to the datasample autocorrelation, where

$\begin{matrix}{{\rho_{\underset{\_}{w}}\lbrack n\rbrack} = \frac{1}{N - {n}}} & (5.33)\end{matrix}$

With these, the redundancy averaged values of Eqn. (5.31) can be writtenasρ_(RA,m) [n]=ρ _(x,m) [n]ρ _(w) [n]  (5.34)

The final values may be averaged over all snapshots,

${\rho_{RA}\lbrack n\rbrack} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}\;{{\rho_{{RA},m}\lbrack n\rbrack}.}}}$The covariance matrix for redundancy averaging is then formed from thesequence asR _(RA)=((ρ_(RA) [r−c]))_(r,c)  (5.35)In terms of the DSR operation,R _(RA)=DSR( R _(SCM))

T _(RA)   (5.36)whereT _(RA)=((ρ _(w) [r−c]))_(r,c)  (5.37)

By way of comparison to the CSS technique discloses herein, the DSRoperation in Eqn. (5.8), which is related to the CSS technique, isdefined as sum, not average. The particular window sampleautocorrelation, ρ _(w) [n]=(N−|n|)⁻¹, used for redundancy averagingresults in an unbiased expected value, E{ρ_(RA)[n]}=R_(x)[n] and E{R_(RA)}=R _(x), if the data is uncorrelated. However, it is also resultsin R _(RA) being an indefinite matrix. This is a familiar result in thecontext of time series analysis, as Eqn. (5.31) is the form of theunbiased estimator for the auto-correlation of a sequence. Whileunbiased, its form in Eqn. (5.34) shows the sequence is a product of twofunctions. The Fourier transform of this product is the estimate of thepower spectral density. It is the convolution

$\begin{matrix}{{{\hat{P}}_{RA}(\psi)} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{C_{\underset{\_}{w}}\left( {\psi - \beta} \right)}{{\hat{P}}_{{\rho\; x},\psi}(\beta)}\ {\mathbb{d}\beta}}}}} & (5.38)\end{matrix}$

where {circumflex over (P)}_(px,ψ)(ψ)=F (ρ_(x)[n]) and C _(w) (ψ)=F (ρ_(w) [n]). For the particular ρw[n] used in redundancy averaging, C _(w)(ψ) is not strictly greater than zero, and as a result portions of{circumflex over (P)}_(RA)(ψ) may become negative valued. This is aninvalid condition for a power spectral density. In its matrix form, thiscondition presents itself by making the covariance R _(RA) indefinite.This is undesirable, and is particular concern for the redundancyaveraging approach.

Covariance Matrix Tapers

Covariance matrix tapers is a technique that provides a measure ofrobustness to sample covariance matrix processing by modifying thesample covariance matrix with a taper matrix, T _(CMT), according toR _(CMT) =R _(SCM) ⊙T _(CMT)  (5.39)

The taper matrix, T_(CMT), is designed specifically to impart nullwidening properties, diagonal loading, or other features. The tapermatrix is positive semi-definite, T _(CMT)≧0, and Hermitian, T _(CMT)=T_(CMT) ^(H). Additionally, the T _(CMT) may be a normalized diagonallyhomogeneous (NDH) matrix, meaning it is a constant down its maindiagonal. Because it is based on the sample covariance matrix, CMT doesnot attempt to address the correlated signal and interference condition.Additional processing, such as spatial smoothing is necessary tomitigate signal cancellation.

Comparison Summary

Thus, the CSS technique disclosed herein provides benefits found in RAand CMT as well as additional benefit. The CSS technique uses the DSRoperation which is beneficial for correlated signal and interference. Asadditional positive attributes, CSS maintains positive definiteness withreasonable restrictions on the choice of window function and resultantpower spectral density estimate. The DSR processing provides additionaldata averaging, or alternatively an increase in effective sample sizeover CMT which uses only the sample covariance matrix.

Non-Ideal Array Manifold Response

Certain CSS methods disclosed herein assume an underlying structurebased upon the observed narrowband space-time process including sums ofphysically propagating plane waves. For an array of idealomnidirectional sensors, the array manifold response takes on a formbased on the complex-exponential, v(k)=((exp[−jk ^(T) p _(n)]))_(n).Real-world sensors and arrays may exhibit perturbations to this idealresponse that alter the form of the encountered covariance from itsassumed structure. For a uniform linear array, the ideal manifoldresponse results in a Toeplitz covariance matrix. With array manifoldresponse error, this Toeplitz structure may no longer hold. This issimilar to the effect encountered when investigating the correlatedsignal and interference scenario. The ML estimate of an unstructuredcovariance is the sample covariance matrix, and the set of algorithmsthat build upon it are most applicable if the underlying problem isunstructured.

Non-ideal array manifold responses and their impact on structuredcovariance beamformer performance are discussed below. A technique tomitigate the impact based on data already available via CSS with MTSEprocessing is described. This is done by estimating the array manifoldresponse corresponding to detectable line components in the spectrum,since these typically dominate the overall performance, andincorporating this discrete set of non-ideal response vectors into thecovariance.

This approach differs from techniques that concentrate on estimating thesteering vector (such as a principal eigenvector of a clutter covariancematrix for radar) or that concentrate on estimates of the actual sensorpositions (such as for towed arrays). Such techniques may use GPS datafor the tow vessel and a water pulley model for the array, andalgorithms to optimize the “focus” or sharpness of the wavenumberspectrum or observations of broadband signals in the environment. As anarrowband processing algorithm operating directly on the snapshot data,the CSS technique disclosed herein may be valuable in concert withlegacy techniques as the array displacement becomes significant, inparticular in light of the similarity between the circular bow arraydeformity and observed array behavior during turning maneuvers.

Types of Non-Ideal Array Manifold Responses

Random Errors

Random errors in array manifold response may be considered zero meanperturbations from the nominal array response values. These may arisefrom non-ideal sensor gain or phase, manufacturing precision, or othercomponent tolerances. Two approaches for describing these effects aredescribed. The first provides a random error term for each physicalquantity related to each sensor, namely, position, amplitude response,and phase response.

Each error term, Δ_(i), is specified as a Gaussian random variable,

${\Delta_{i}\overset{d}{\longrightarrow}{N\left( {0,\sigma_{i}^{2}} \right)}}.$The “actual” value of a quantity, x, may be denoted using a subscript a(e.g., x_(a)). The effective position of the n^(th) sensor element withnominal position p_(n) isp _(n,a) =p _(n)+[Δ_(x),Δ_(y),Δ_(z)]^(T)  (6.1)

and the overall array manifold response isv _(a)=(([1+Δ_(A)]exp[jΔ _(θ)]exp[−j( k ^(T) p _(n,a))]))_(n)  (6.2)

This provides complete specification of the potential random errors,(Δ_(x), Δ_(y), Δ_(z), Δ_(A), Δ_(θ)), and is useful for analysis ofperformance impacts in relation to each of the individual quantities. Asimpler model reduces the random error contribution to a multiplicativeeffect on the amplitude and phase only. Define the vector

$\begin{matrix}{{\underset{\_}{h} = {{c_{h}\underset{\_}{1}} + \underset{\_}{g}}},{\underset{\_}{g}\overset{d}{\longrightarrow}{{CN}_{N}\left( {0,{\sigma_{g}^{2}\underset{\_}{I}}} \right)}}} & (6.3)\end{matrix}$

The non-ideal array manifold response isv _(a)= v

h =c_(h) v + v _(g)  (6.4)

The constant c_(h) is chosen such thatv _(a) ^(H) v _(a) =v ^(H) v   (6.5)

The ratio σ_(g) ²/c_(h) ² is a single metric that provides a measure ofthe difference between the ideal and actual array manifold responses.Expressed in dB as 10 log₁₀(σ_(g) ²/c_(h) ²), this value gives someindication of “how far down” the perturbation components are from thenominal response. A simpler model may be useful for varying arrayresponse as a function of angle of arrival, while the more explicitmodel may be more appropriate for analyzing operation below the designfrequency where there is a non-zero virtual region and isotropic noisecomponent.

Deterministic Errors

In addition to random errors, arrays may experience more deterministictypes of array manifold response perturbation. For example, deformationof the array may cause a non-zero mean, non-random disturbance in thepositions of the elements. With underwater acoustic arrays, suchdeformation can occur due to hydrodynamics for towed arrays in motion.Two types of positional errors for linear arrays, circular and partialcircular bows, are considered below.

Impact to Structured Covariance Matrix ABF Performance

The technique of computing covariance based on estimates of thewavenumber spectrum may be influenced by factors that would impact theability to accurately estimate that spectrum, or how accurately thatspectrum (based on complex exponentials) represents the true underlyingsituation. Using the simple model described above, the actual arraymanifold response vector is v _(a)=c_(h) v+v _(g). The true arraymanifold response, v, is scaled and the additional bias term is v_(g)=v⊙g. Because v(k)=((e^(−jk) ^(T) ^(p) ^(n) ))_(n):

$\begin{matrix}{\underset{\_}{g},{{\underset{\_}{v}}_{g}\overset{d}{\longrightarrow}{{CN}_{N}\left( {0,{\sigma_{g}^{2}\underset{\_}{I}}} \right)}}} & (6.6)\end{matrix}$

so that the bias in any given realization is a complex Gaussian randomvector. Consider a single interferer in uncorrelated white noise.x _(m) =v _(a) a(m)+ n _(m)  (6.7)

The covariance for this case for a given instance of v _(g) isE{x _(m) x _(m) ^(H) |vg}=σ ²(|c _(h)|² vv ^(H) +c _(h) vv _(g) ^(H) +c_(h) *v _(g) v ^(H) +v _(g) v _(g) ^(H))+ R _(n)  (6.8)

In any given instance of this scenario, the error is unknown butnon-random, and will produce a particular wavenumber spectrum. Thecovariance over the ensemble of error vectors, v _(g), isE{x _(m) x _(m) ^(H)}=σ² |c _(h)|² vv ^(H) +R _(n)+σ²σ_(g) ² I   (6.9)

The ensemble power spectrum for this case appears as the original withthe line component scaled | c_(h)|², with an elevated noise floorcorresponding to the σ²σ_(g) ² I term. In a given realization though, itis a single vector, v _(g), causing the elevated sidelobes and not anensemble of plane waves from all directions. As the ratio σ_(g)/c_(h) ²increases there is larger and larger array manifold error, and the risein the perceived noise floor is clear.

The normalized SINR loss performance of structured covariance techniquesmay degrade due to random errors and may be proportional to thedifference between the sidelobe levels due to the errors and the truenoise floor. The expected increase in normalized SINR loss is equal tothe rise in the noise floor.Δξ_(dB)≈max(0,INR_(dB)+(σ_(g) ² /c _(h) ²)dB)  (6.10)

Deterministic array manifold response error also causes an impact to theestimated wavenumber spectrum. This impact is a function of both thepositional displacement and the angle of arrival of the discrete pointsource. Due to the symmetry of the array bend, the sidelobe structureresulting from the circular bow is expected to be symmetric inwavenumber space. The circular bow distortion causes the single discretesource to appear as a symmetric, spatially spread interference.

Mitigation Techniques

Non-ideal array manifold response may be expressed in a snapshot modelby adding an error vector (or calibration error), u _(k), to the idealresponse vector, v _(k), in building the snapshots.

$\begin{matrix}{{\underset{\_}{x}}_{m} = {{\sum\limits_{k = 1}^{K}\;{\left( {{\underset{\_}{v}}_{k} + {\underset{\_}{u}}_{k}} \right){a_{k}(m)}}} + {\underset{\_}{n}}_{m}}} & (6.11)\end{matrix}$

The impact of combinations of array manifold response errors and thebackground noise, n _(m),are not considered since the dominant source ofperformance degradation is assumed to be caused by the large INR linecomponents. Because the error vectors are additive to the ideal responsevector, the snapshots include a set of terms corresponding to the idealarray manifold response and the additional error terms.

$\begin{matrix}{{\underset{\_}{x}}_{m} = {\left\lbrack {{\sum\limits_{k = 1}^{K}\;{{\underset{\_}{v}}_{k}{a_{k}(m)}}} + {\underset{\_}{n}}_{m}} \right\rbrack + {\sum\limits_{k = 1}^{K}\;{{\underset{\_}{u}}_{k}{a_{k}(m)}}}}} & (6.12)\end{matrix}$

The ideal array manifold response drives the structure in the CSSalgorithms. The case considered here is where the array manifoldresponse errors exist, but are small in magnitude compared to the truearray manifold response. This situation is referred to herein as“partially structured”, to reflect that the ideal array manifoldresponse is still somewhat discernible in the data. The MTSE harmonicanalysis still provides the estimates of the observable plane wavecomponents, provided they maintain a sufficient SNR over the apparentraised noise floor such that line component detection is possible.

The random error case may be considered as a starting point to developan algorithm to deal with random array manifold perturbation. In aparticular embodiment, an optimal algorithm to deal with deterministicarray manifold errors may estimate the error model parameters, H for acircular bow, and H and L2 for a partial circular bow, or parabolicparameters in, and may use that information to assist in estimating thearray manifold errors. The discussion below further describes aparticular technique developed for random errors, and considers to whatextent deterministic errors may be processed effectively for two typesof non-random errors that are considered.

MMSE Unbiased Linear Estimate of Calibration Errors

A minimum mean squared error, linear unbiased estimate of calibrationerrors is explained below. For ease of explanation, there are assumed tobe K line components, and M≧K snapshots observed. A snapshot modelcontaining non-ideal array manifold responses, v _(k,a)=v _(k)+u _(k),provides a starting point. The collection of snapshots, m=1, . . . , Mcan be arranged as a matrixX=[x ₁ ,x ₂ , . . . ,x _(M)]=( V+U ) A+N   (6.13)

where the individual matrices similarly include the original vectors,V=[v ₁, v ₂, . . . , v _(K)], U=[u ₁, u ₂, . . . , u _(K)], A=[a ₁, a ₂,. . . , a _(M)], and N=[n ₁, n ₂, . . . , n _(M)]. Here, a prioriknowledge of the ideal array manifold responses, V, and signalamplitudes, A, may be assumed. In practice, the estimates provided bythe MTSE harmonic analysis processing may be used for these values. Thematrix A may be represented in an additional way, with each rowrepresenting the amplitude time series for the corresponding pointsource signal.α _(k) ^(T) =[a _(k)(1),a _(k)(2), . . . ,a _(k)(M)]  (6.14)A=[α ₁,α ₂, . . . ,α _(K)]^(T)  (6.15)

The array calibration errors U may be estimated given V and A. Aninitial step may be performed to subtract out the known line components.Y=X−VA=UA+N   (6.16)

Additional calculations may be performed subject to the followingassumptions:

The array calibration errors, u _(k), are non-random but unknown and aredifferent for each source, k=1 . . . K.

The noise terms are zero mean and independent of the array manifoldresponses and signal amplitudes, with covariance, E {n _(m) n _(m)^(H)}=R _(n).

The snapshots for different sample indices, m, are independent.

rank(A)=K, which is discussed further below. This rank deficientsituation may occur when two signals are perfectly correlated with oneanother (the magnitude of the correlation coefficient is unity), or whenthe number of interferers exceeds the number of snapshots, K>M.

In a particular embodiment, a linear processor can be used to estimatethe array manifold response error vectors from the observations. Foreach error vector (or for all û _(k) simultaneously){circumflex over (u)} _(k) =Yw _(k),{circumflex over (U)}= YW   (6.17)

where the matrices are W=[w ₁, w ₂, . . . , w _(K)] and Û=[û ₁, û ₂, . .. , û _(K)]. The processor may be unbiased, so that E {û _(k)}=u _(k),E={Û}U. Expanding out the expectationE{{circumflex over (u)} _(k) }=E{Yw _(k) }=UAw _(k)  (6.18)

This implies that to be unbiased UAw _(k)=u _(k), orAw _(k) =e _(k), or AW=I   (6.19)

where e _(k) is the elementary vector which is all zeros with a 1 in thek^(th) position. This places K constraints on each weight vector, w_(k). Written out explicitly using a row definition of Aα _(j) ^(T) w _(k) =δ[j−k]  (6.20)

where δ[j−k] is a Kronecker delta. In a particular embodiment, theestimation error variance may be minimized at the output for eachfilter, w _(k), given byσ _(w,k) ² =E[( u _(k)−{circumflex over (u)} _(k))^(H)( u_(k)−{circumflex over (u)} _(k))]  (6.21)

Using Eqn. (6.17) in Eqn. (6.21), the expression for the variancereduces toσ _(w,k) ² =w _(k) ^(H) E[N ^(H) N]w _(k)  (6.22)

To find E {N ^(H) N} recall that N=[n ₁, n ₂, . . . , n _(M)], where theindividual n _(m) are i.i.d. with covariance, R _(n).

$\begin{matrix}\begin{matrix}{{E\left\{ {{\underset{\_}{N}}^{H}\underset{\_}{N}} \right\}} = {E\left\{ \begin{pmatrix}{{\underset{\_}{n}}_{1}^{H}{\underset{\_}{n}}_{1}} & {{\underset{\_}{n}}_{1}^{H}{\underset{\_}{n}}_{2}} & \ldots & {{\underset{\_}{n}}_{1}^{H}{\underset{\_}{n}}_{M}} \\{{\underset{\_}{n}}_{2}^{H}{\underset{\_}{n}}_{1}} & {{\underset{\_}{n}}_{2}^{H}{\underset{\_}{n}}_{2}} & \; & \; \\\vdots & \; & \ddots & \; \\{{\underset{\_}{n}}_{M}^{H}{\underset{\_}{n}}_{1}} & \; & \; & {{\underset{\_}{n}}_{M}^{H}{\underset{\_}{n}}_{M}}\end{pmatrix} \right\}}} \\{= {{{tr}\left( {\underset{\_}{R}}_{n} \right)}\underset{\_}{I}}}\end{matrix} & (6.23)\end{matrix}$

Inserting Eqn. (6.23) into Eqn. (6.22) gives the noise output power.σ _(w,k) ² =tr( R _(n)) w _(k) ^(H) w _(k)  (6.24)

The constrained optimization problem may be set up to find theindividual weights, w _(k).

$\begin{matrix}{{\begin{matrix}{\arg\;\min} \\{\underset{\_}{w}}_{k}\end{matrix}{{tr}\left( {\underset{\_}{R}}_{n} \right)}{\underset{\_}{w}}_{k}^{H}{\underset{\_}{w}}_{k}\mspace{14mu}{subject}\mspace{14mu}{to}\mspace{14mu}{\underset{\_}{\alpha}}_{j}^{T}{\underset{\_}{w}}_{k}} = \delta_{jk}} & (6.25)\end{matrix}$

The following cost function may be minimized using the method ofLagrange multipliers.

$\begin{matrix}{{J\left( {\underset{\_}{w}}_{k} \right)} = {{{{tr}\left( {\underset{\_}{R}}_{n} \right)}{\underset{\_}{w}}_{k}^{H}{\underset{\_}{w}}_{k}} + {\sum\limits_{l = 1}^{K}\;{\lambda_{kl}\left( {{{\underset{\_}{\alpha}}_{l}^{T}{\underset{\_}{w}}_{k}} - \delta_{lk}} \right)}} + {\sum\limits_{l = 1}^{K}\;{\lambda_{kl}^{*}\left( {{{\underset{\_}{w}}_{k}^{H}{\underset{\_}{\alpha}}_{l}^{*}} - \delta_{lk}} \right)}}}} & (6.26)\end{matrix}$

Taking the gradient w.r.t. w _(k) and setting to 0

$\begin{matrix}{{\frac{\partial}{\partial{\underset{\_}{w}}_{k}}{J\left( {\underset{\_}{w}}_{k} \right)}} = \left. 0 \right|_{{\underset{\_}{w}}_{k} = {\underset{\_}{w}}_{k,{opt}}}} & (6.27)\end{matrix}$

yields the following expression

$\begin{matrix}{{\underset{\_}{w}}_{k}^{H} = {{- \frac{1}{t\;{r\left( {\underset{\_}{R}}_{n} \right)}}}{\sum\limits_{l = 1}^{K}{\lambda_{k\; l}{\underset{\_}{\alpha}}_{l}^{T}}}}} & (6.28)\end{matrix}$

Defining the vector λ _(k) and matrix Λ respectively,λ _(k)=[λ_(k,1),λ_(k,2), . . . ,λ_(k,K)],Λ=[λ ₁ ^(T),λ ₂ ^(T), . . . ,λ_(K) ^(T)]^(T)  (6.29)

then Eqn. (6.28) can be written as

$\begin{matrix}{{\underset{\_}{w}}_{k}^{H} = {{- \frac{1}{t\;{r\left( {\underset{\_}{R}}_{n} \right)}}}{\underset{\_}{\lambda}}_{k}\underset{\_}{A}}} & (6.30)\end{matrix}$

The entire solution can be expressed for all k simultaneously as

$\begin{matrix}{{\underset{\_}{W}}^{H} = {{- \frac{1}{t\;{r\left( {\underset{\_}{R}}_{n} \right)}}}\underset{\_}{\Lambda}\underset{\_}{A}}} & (6.31)\end{matrix}$

Inserting Eqn. (6.31) into the constraint Eqn. (6.19), and solving forthe constraint matrix, ΛΛ=−tr(R _(n))(AA ^(H))⁻¹  (6.32)

The earlier assumption that rank(A)=K implies that (AA ^(H))⁻¹ exists.Inserting Eqn. (6.32) back into Eqn. (6.31) produces the overallsolution for the weights.W=A ^(H)( AA ^(H))⁻¹  (6.33)

which is the Moore-Penrose pseudo-inverse of A. Using Eqn. (6.33) inEqn. (6.17), the minimum variance unbiased linear estimate of the nonideal array manifold response error vectors is{circumflex over (U)}= YW=YA ^(H)( AA ^(H))⁻¹  (6.34)

In situations where M<K, then rank(A)<K and the inverse (AA ^(H))⁻¹ doesnot exist. In this situation, the processing may be restricted toestimate a subset of the array manifold error vectors, û _(k),corresponding to the M largest sources. In practice, there may bepractical issues with the conditioning of the matrix AA ^(H) and a valueless than M may work better numerically. Performance may degrade inthese scenarios compared to an optimal adaptive beamformer because thereis not enough snapshot information to estimate all the random errorvectors from the data.

Maximum Likelihood Estimate of Calibration Errors

In the discussion above, the minimum mean square error linear unbiasedestimate of the array manifold response error vectors was determined.This used knowledge of the number of point source signals, K, the signalamplitudes, A, and the ideal array manifold response vectors, V. Thesnapshots were assumed to be independent. The solution did not requiredefinition of the statistics of the noise term, n _(m), other than tostate that the noise was independent of the other quantities in thedata. In the discussion below, the maximum likelihood estimation of thearray manifold response error vectors is derived. The same assumptionsare used as were used above. Additionally, the noise terms are specifiedas complex Gaussian,

$\underset{\_}{n}\overset{d}{\rightarrow}{C\;{{N_{N}\left( {\underset{\_}{0},{\underset{\_}{R}}_{n}} \right)}.}}$The individual snapshotsx _(m)=( V+U ) a _(m) +n _(m)  (6.35)

are complex Gaussian random vectors, with non-zero mean, m _(m)=[V+U]a_(m)

$\begin{matrix}{{\underset{\_}{x}}_{m}\overset{d}{\rightarrow}{C\;{N_{N}\left( {{\underset{\_}{m}}_{m},{\underset{\_}{R}}_{n}} \right)}}} & (6.36)\end{matrix}$

Grouping all snapshots, X=[x ₁, x ₂, . . . , x _(M)] and amplitudevalues, A=[a ₁, a ₂, . . . , a _(M)], the complex Gaussian random matrixX has non-zero mean M=(V+U) A. With the snapshots being independent, theprobability density for X is

$\begin{matrix}{{f_{\underset{\_}{x}}\left( \underset{\_}{X} \right)} = {\pi^{{- \; N}\; M}{\det\left( {\underset{\_}{R}}_{n} \right)}^{- M}{\exp\left( {- {\sum\limits_{m = 1}^{M}{\left( {{\underset{\_}{x}}_{m} - {\underset{\_}{m}}_{m}} \right)^{H}{{\underset{\_}{R}}_{n}^{- 1}\left( {{\underset{\_}{x}}_{m} - {\underset{\_}{m}}_{m}} \right)}}}} \right)}}} & (6.37)\end{matrix}$

Eqn. (6.37) can be expressed using the trace operator, tr( ).f _(x) ( X )=π^(−NM) det( R _(n))^(−M)exp(−tr[( X - M )^(H) R _(n) ⁻¹(X - M )])  (6.38)

For the available snapshot data, the likelihood function for a givenestimate of the array manifold errors isf _(X|Û) ( X |{circumflex over (U)})π^(−NM) det( R _(n))^(−M)exp(−tr[(X - M ({circumflex over ( U )}))^(H) R _(n) ⁻¹( X - M ({circumflex over(U)}))])  (6.39)

The maximum likelihood estimate of the array manifold error vectors,Û_(ML), is a matrix that maximizes f _(X|Û) (X|Û).

$\begin{matrix}{{\underset{\_}{\hat{U}}}_{ML} = {\begin{matrix}{\arg\;\max} \\\underset{\_}{\hat{U}}\end{matrix}{f_{X|\hat{\underset{\_}{U}}}\left( \underset{\_}{X} \middle| \hat{\underset{\_}{U}} \right)}}} & (6.40)\end{matrix}$

Beginning with the log-likelihood function,

_(X|Û) (X|Û)=ln f _(X|Û) (X|Û).

X _ | U ^ _ ⁢ ( X _ | U ^ _ ) = C 1 + C 2 - tr ⁡ [ ( X _ - M _ ⁡ ( U ^ _ )) H ⁢ R _ n - 1 ⁡ ( X _ - M _ ⁡ ( U ^ _ ) ) ] ( 6.41 )

where the constants C₁, C₂ are not functions of Û. Because it is amaximum of

_(X|Û) (X|Û).

∂ ∂ U ^ _ ⁢ X _ | U ^ _ ⁢ ( X _ | U ^ _ ) ⁢ | U ^ _ = U ^ _ M ⁢ ⁢ L = 0 (6.42 )

Expanding out terms in Eqn. (6.41), and applying the partial derivativethrough the tr ( ) function{circumflex over (U)} _(ML) =XA ^(H)( AA ^(H))⁻¹ −V   (6.43)

Alternatively, starting withY=X−VA   (6.44)

and applying the same procedure{circumflex over (U)} _(ML) =YA ^(H)( AA ^(H))⁻¹  (6.45)

Eqn. (6.45) matches the result Eqn. (6.34). It is interesting to note isthat the covariance of the noise term n _(m) does not have an impact onthe final solution. This is particularly useful, since this covarianceis not known at this point in the processing. Also, it establishes thatproblems should not be encountered if the underlying noise is non-white,due to either spatially spread interference or the combination ofvisible/virtual space when operating below the design frequency.

Where a non-white noise may still cause issues is in situations wherethe dynamic range of the noise, spectrally, is very large. This could bethe case when operating below design frequency where the isotropic noisecomponent is much larger than the sensor noise component. The apparentincrease in noise floor due to the effects of array manifold responseerrors accumulated by the spatially spread isotropic noise process maymask the true sensor noise components.

Procedure

The discussion below describes a procedure to use the results fromderivation of Eqns. (6.45) and (6.34) to supplement the CSS with MTSEapproach with additional non-ideal array manifold response error vectorinformation. The dominant sources of normalized SINR loss are the errorvectors associated with the high INR point sources in the data. The termpartially structured is used to indicate that some observation of theunderlying structure remains observable. Fundamentally, being partiallystructured implies that the high INR line components in the spectrum canstill be detected via the MTSE harmonic analysis process. In general,this indicates some amount of positive SNR of the line component abovethe raised noise floor caused by the array manifold errors.

The following procedure may be used to incorporate non-ideal arraymanifold response error vector information into an overall covariancematrix estimate.

1. Use MTSE harmonic analysis to detect and estimate the number of linecomponents, K, and their parameters, v _(k), and a_(k)(m).

2. Use Eqns. (6.44) and (6.45) to estimate array calibration errorvectors, u _(k), for these line components.

3. Subtract the discrete components from the snapshot data to form theresidual, X _(res)=X−(V+U)A.

4. Optionally, items 1, 2 and 3 may be iterated.

5. Estimate the final residual “continuous” background spectrum,{circumflex over (P)}_(res)(·), and form R _(MTSE,res).

6. Form an overall estimate of the covariance matrix

$\begin{matrix}{{\underset{\_}{R}}_{TOTAL} = {{\sum\limits_{k = 1}^{K}{\left( {{\underset{\_}{v}}_{k} + {\underset{\_}{u}}_{k}} \right)\left( {{\underset{\_}{v}}_{k} + {\underset{\_}{u}}_{k}} \right)^{H}\frac{1}{M}{\sum\limits_{m = 1}^{M}{{a_{k}(m)}}^{2}}}} + {\underset{\_}{R}}_{{MTSE},{res}}}} & (6.46)\end{matrix}$

Comparison-MTSE and MVDR Spectra

The error vector processing is a technique for adapting the CSS withMTSE approach to handle non-ideal array manifold response, where theresponse retains sufficient “likeness” to the ideal that the parametersof the large discrete components in the spectrum can be estimated. Usingthis information, the error vector(s) can be estimated and used in theformulation of the final estimate of covariance. Simulations that havebeen performed indicate that this processing is effective and that theoverall normalized SINR loss is equivalent to reduced rank processing ofthe data using MWF. Array manifold error processing may not be areplacement for unconstrained covariance and reduced rank techniques inall circumstance, but it offers a way to extend the basic CSS with MTSEprocessing based on available data products to address the issue ofnon-ideal response.

An additional byproduct of this processing is a final MTSE wavenumberspectrum with the influence of the error vectors for the large discretesremoved. One can qualitatively compare the utility in this finalestimate of the spectrum to MVDR spectra computed using the ensemblecovariance matrix.

The techniques described above may be used to estimate the error vectorsassociated with the strongest INR line components in the data using thedetection and estimation parameters available from harmonic analysis.Through simulations, this technique was seen to improve performance inthe face of non-ideal array manifold response, and normalized SINR losswas seen to be comparable to reduced rank techniques provided theharmonic analysis could reliably detect the line components. As abeneficial byproduct, MTSE could be made to produce an estimate of thewavenumber spectrum as if the array manifold had been restored to ideal.

Extensions for Arbitrary Geometry

Arrays with regular element spacing lend themselves to a convenientsimultaneous estimation of visible space components of the space-timeprocess and virtual space sensor noise using classical spectralestimation techniques and efficient FFT computation. For arbitrary arraygeometries, proper estimation of the virtual space sensor noise mayrequire more effort. The following discussion explains an approach fordoing this based upon analyzing the covariance for isotropic noise.Methods to extend the disclosed CSS with MTSE techniques for applicationto arbitrary arrays are also explained below.

Sensor Noise with Arbitrary Array Geometry

The covariance matrix of interest, R _(x), includes two subspacescorresponding to visible and virtual regions.< R _(x) >=<R _(vs) >+<R _(vr)>  (7.1)

Assuming that the subspaces, <R _(vs)> and <R _(vr)> are approximatelyorthogonal, although due to the array geometry the transition betweenthe two may not be so sharp. The visible region subspace isapproximately the subspace defined by the covariance associated with 3Disotropic noise.< R _(vs) >≈<R _(iso)>  (7.2)

This is convenient because the 3D isotropic noise is specified in termsof angle of arrival to the array. This has an intuitive physicalinterpretation regardless of array geometry. Additionally, from Eqn.(3.31), the contribution of the sensor noise component appeared as a 3Disotropic noise component when restricting attention to the visibleregion. This may be useful later when considering the positivedefiniteness of the estimated covariance.

Because of the spatial stationarity of the 3D isotropic noise, thecovariance is a function of the difference in position, Δp. Thisdifference may be expressed in Cartesian coordinates {Δp_(x), Δp_(y),Δp_(z)} or spherical coordinates {s, γ, ζ}.

$\begin{matrix}{{\Delta\underset{\_}{p}} = {\begin{pmatrix}{\Delta\; p_{x}} \\{\Delta\; p_{y}} \\{\Delta\; p_{z}}\end{pmatrix} = \begin{pmatrix}{s\mspace{14mu}\sin\mspace{14mu}\gamma\mspace{14mu}\cos\mspace{14mu}\zeta} \\{s\mspace{14mu}\sin\mspace{14mu}\gamma\mspace{14mu}\sin\mspace{14mu}\zeta} \\{s\mspace{14mu}\cos\mspace{14mu}\gamma}\end{pmatrix}}} & (7.3)\end{matrix}$

The covariance between two omnidirectional sensors in 3D isotropic noisehas a form

$\begin{matrix}{{R_{iso}\left( {s,\gamma,\zeta} \right)} = {{R_{iso}(s)} = \frac{\sin\left( {\omega\;{s/c}} \right)}{\omega\;{s/c}}}} & (7.4)\end{matrix}$

The covariance matrix for an array of sensors is populated by values ofthis function where the relative position is Δp=p _(r)−p _(c).R _(iso)=(( R _(iso)(Δ p )))_(r,c)=((R _(iso)(s)))_(r,c)  (7.5)

Eqn. (7.5) may be used to compute R _(iso), and to perform aneigendecomposition of the matrix. The decomposition is ordered accordingto decreasing eigenvalues.

$\begin{matrix}{{\underset{\_}{R}}_{iso} = {{{\underset{\_}{Q}}_{iso}{\underset{\_}{\Lambda}}_{iso}{\underset{\_}{Q}}_{iso}^{H}} = {\sum\limits_{n = 0}^{N - 1}{\lambda_{n}{\underset{\_}{q}}_{n}{\underset{\_}{q}}_{n}^{H}}}}} & (7.6)\end{matrix}$

The eigenvalues for this problem, λ_(n), are a measure of theconcentration of each eigenvector within the visible region subspace.The eigenvectors may be divided into two sets, those mostly contained inthe visible region subspace and those mostly contained in the virtualregion subspace. The boundary may be apparent from inspection of theeigenvalues.

Indicating the number of eigenvectors determined to be in the visiblespace as N_(vs), show the explicit make up of R _(iso).

$\begin{matrix}{{\underset{\_}{R}}_{iso} = {{\sum\limits_{n = 0}^{N_{v\; s} - 1}{\lambda_{n}{\underset{\_}{q}}_{n}{\underset{\_}{q}}_{n}^{H}}} + {\sum\limits_{n = N_{v\; s}}^{N - 1}{\lambda_{n}{\underset{\_}{q}}_{n}{\underset{\_}{q}}_{n}^{H}}}}} & (7.7)\end{matrix}$

The visible region subspace is approximately the subspace spanned by thefirst N_(vs) eigenvectors.< R _(vs) >≈<q ₀ ,q ₁ , . . . ,q _(N) _(vs) ⁻¹)>  (7.8)

The eigenvectors, Q _(iso), make a complete orthonormal set for<C^(N×N)>. The covariance matrix R _(x) is full rank, so from<R_(x)><C^(N×N)>=<Q _(iso)>. Thus,< R _(vr) >≈<q _(N) _(vs) ,q _(N) _(vs) ₊₁ , . . . ,q _(N-1)>  (7.9)

Grouping the appropriate eigenvectors together for the visible andvirtual region subspaces,< Q _(vs) =[q ₀ ,q ₁ , . . . ,q _(Nvs-1) ],Q _(vr) =[q _(Nvs) ,q_(Nvs+1) , . . . ,q _(N-1)],  (7.9)

projection matrices for those subspaces may be formedP _(vs) =Q _(vs) Q _(vs) ^(H) ,P _(vr) =Q _(vr) Q _(vr) ^(H) =I - Q_(vs) Q _(vs) ^(H)  (7.11)

Let N_(vr)=N−N_(vs). The sensor noise power component within the virtualregion subspace can be estimated from the available snapshots as

$\begin{matrix}{{\hat{\sigma}}_{w}^{2} = {{\frac{1}{M} \cdot \frac{1}{N_{v\; r}}}{\sum\limits_{m = 1}^{M}{{\underset{\_}{x}}_{m}^{H}{\underset{\_}{P}}_{v\; r}{\underset{\_}{x}}_{m}}}}} & (7.12)\end{matrix}$

From the form of the projection matrix, P _(vr), the sensor noise powerestimate in Eqn. (7.12) is equivalent to

$\begin{matrix}{{\hat{\sigma}}_{w}^{2} = {{\frac{1}{M} \cdot \frac{1}{N_{vr}}}{\sum\limits_{m = 1}^{M}{\sum\limits_{n = N_{vs}}^{N - 1}{{{\underset{\_}{x}}_{m}^{H}{\underset{\_}{q}}_{n}}}^{2}}}}} & (7.13)\end{matrix}$

Eqn. (7.13) shows that {circumflex over (σ)}_(w) ² is the average poweracross each of the orthonormal basis vectors, q _(n), in the subspace.The overall estimate for the covariance R _(x) then uses this estimateas{circumflex over (R)} _(x)={circumflex over (R)} _(vs)+{circumflex over(σ)}_(w) ² P _(vr)  (7.14)

Use of the projection matrix P _(vr) in Eqn. (7.14) avoids doublecounting the sensor noise component measured in the visible regionsubspace and contained in {circumflex over (R)}_(vs). If the sensornoise is significantly below the continuous background noise componentof the space-time process throughout the visible region, one could usethe simpler expression{circumflex over (R)} _(x)={circumflex over (R)} _(vs)+{circumflex over(σ)}_(w,vr) ² I   (7.15)

at the expense of double counting the sensor noise in the visibleregion. In either case some representation of the sensor noise in thevirtual region may be used.

Positive Definiteness

The estimated covariance must be positive definite to be meaningful forarray processing. For regularly spaced arrays, the estimated covarianceis positive definite when the estimated frequency-wavenumber spectrum,{circumflex over (P)}(k)>0, and the estimated sensor noise, {circumflexover (σ)}_(w) ²>0, as explained above. The estimated covariance is alsopositive definite for arbitrary geometry arrays, where the covariance isestimated as the sum of the visible space covariance, {circumflex over(R)} _(vs), and the sense noise in the virtual space.

Method 1

Eqn. (7.14) provides an accurate method for accounting for the sensornoise in the covariance matrix.{circumflex over (R)} _(x)={circumflex over (R)} _(vs)+σ_(w) ² P_(vr)  (7.16)

The estimate for the covariance corresponding to the visible region ofthe array is

$\begin{matrix}{{\underset{\_}{\hat{R}}}_{vs} = {\left( {2\pi} \right)^{- C}{\int{\ldots{\int_{vs}{{{\hat{P}}_{vs}\left( \underset{\_}{k} \right)}{\underset{\_}{v}\left( \underset{\_}{k} \right)}{{\underset{\_}{v}}^{H}\left( \underset{\_}{k} \right)}{\mathbb{d}\underset{\_}{k}}}}}}}} & (7.17)\end{matrix}$

The estimate of the wavenumber spectrum in the visible region may beassumed to be greater than zero, {circumflex over (P)}_(vs)(k)>0, whichenables expressing it as including an estimate of the wavenumberspectrum of the observed process, {circumflex over (P)}_(f)(k)≧0, and anestimate of the sensor noise seen in the visible region, {circumflexover (σ)}_(w,vs) ²>0.{circumflex over (P)} _(vs)( k )={circumflex over (P)} _(f)( k)+{circumflex over (σ)}_(w,vs) ²  (7.18)

Using Eqns. (7.18) and (7.17) in the quadratic expression for positivedefinitenessx ^(H) {circumflex over (R)} _(x) x=x ^(H)((2π)^(−C)∫ . . . ∫_(vs)[{circumflex over (P)} _(f)( k )+{circumflex over (σ)}_(w,vs) ² ]v ( k )v ^(H)( k )dk+{circumflex over (σ)} _(w) ² P _(vr)) x   (7.19)

Carrying out the integration, this simplifies tox ^(H) {circumflex over (R)} _(x) x=x ^(H) {circumflex over (R)} _(f)x+x ^(H)({circumflex over (σ)}_(w,vs) ² R _(iso)+σ_(w) ² P _(vr)) x  (7.20)

The covariance estimate for the space-time process is positivesemi-definite, {circumflex over (R)}_(f)≧0. This follows fromx ^(H) {circumflex over (R)} _(f) x =(2π)^(−C)∫ . . . ∫_(vs) {circumflexover (P)} _(f)( k )| x ^(H) v ( k )|² dk   (7.21)

where {circumflex over (P)}_(f)(k)≧0 and |x ^(H) v(k)|²≧0. For the termsrelating to the noise estimates, R _(iso) may be replaced with itseigendecomposition.

$\begin{matrix}{{{{\hat{\sigma}}_{w,{vs}}^{2}{\underset{\_}{R}}_{iso}} + {{\hat{\sigma}}_{w}^{2}{\underset{\_}{P}}_{vr}}} = {{{\hat{\sigma}}_{w,{vs}}^{2}\left\lbrack {\sum\limits_{n = 0}^{N - 1}{\lambda_{n}{\underset{\_}{q}}_{n}{\underset{\_}{q}}_{n}^{H}}} \right\rbrack} + {\sigma_{w}^{2}{\underset{\_}{P}}_{v\; r}}}} & (7.22)\end{matrix}$

The matrix P _(vr) may also be defined in terms of theeigendecomposition of R _(iso)

$\begin{matrix}{{{{\hat{\sigma}}_{w,{vs}}^{2}{\underset{\_}{R}}_{iso}} + {{\hat{\sigma}}_{w}^{2}{\underset{\_}{P}}_{vr}}} = {{\sum\limits_{n = 0}^{N_{vs} - 1}{{\hat{\sigma}}_{w,{vs}}^{2}\lambda_{n}{\underset{\_}{q}}_{n}{\underset{\_}{q}}_{n}^{H}}} + {\sum\limits_{n = N_{vs}}^{N - 1}{\left( {{{\hat{\sigma}}_{w,{vs}}^{2}\lambda_{n}} + {\hat{\sigma}}_{w}^{2}} \right){\underset{\_}{q}}_{n}{\underset{\_}{q}}_{n}^{H}}}}} & (7.23)\end{matrix}$

Define the combined eigenvalues, λ_(n)′, as

$\begin{matrix}{\lambda_{n}^{\prime} = \left\{ \begin{matrix}{{\hat{\sigma}}_{w,{vs}}^{2}\lambda_{n}} & {0 \leq n \leq N_{vs}} \\{{{\hat{\sigma}}_{w,{vs}}^{2}\lambda_{n}} + {\hat{\sigma}}_{w}^{2}} & {N_{vs} \leq n < N}\end{matrix} \right.} & (7.24)\end{matrix}$

By the selection methods outlined above, the eigenvalues λ_(n) for0≦n<N_(vs) are greater than zero, such that all λ_(n)′ are real-valuedand greater than zero. Defining Λ′=diag (λ₀′, λ₁′, . . . , λ′_(N-1))x ^(H)({circumflex over (σ)}_(w,vs) ² R _(iso)+σ_(w) ² P _(vr)) x=x ^(H)Q _(iso) Λ′Q _(iso) ^(H) x   (7.25)

The matrix Q_(iso) is unitary, therefore |Q _(iso) ^(H) x|²=|x|²≠0 for|x|²≠0. With positive, real-valued entries on the diagonal of Λ′x ^(H)({circumflex over (σ)}_(w,vs) ² R _(iso)+{circumflex over (σ)}_(w)² P _(vr)) x>0  (7.26)

The overall matrix {circumflex over (R)} _(x) is the sum of a positivesemidefinite matrix, {circumflex over (R)} _(f), and a positive definitematrix, x ^(H) ({circumflex over (σ)}_(w,vs) ² R _(iso)+{circumflex over(σ)}_(w) ² P _(vr))x, and therefore is positive definite, {circumflexover (R)} _(x)>0.

Method 2

Another method for an arbitrary array geometry is{circumflex over (R)} _(x)={circumflex over (R)} _(vs)+{circumflex over(σ)}_(w,vs) ² I   (7.27)

As explained above, {circumflex over (R)} _(vs)=({circumflex over (R)}_(f)+{circumflex over (σ)}_(w,vs) ² R _(iso))≧0. The matrix I>0 and so{circumflex over (R)} _(x)>0.

Covariance from Spatial Spectra (CSS) with MTSE

Design of Multi-Tapers

For arbitrary array geometry, the design of the multiple tapers forspectral estimation may be done by dividing the visible region into asearch grid, where each grid point defines a region of analysis. Forconvenience, the angle domain, (θ, φ), is used herein, althoughoperation may alternately be specified in wavenumber. In general form,the search grid covers a sphere, 0≦θ≦π, 0≦φ≦2π, although the arraycharacteristics may be exploited to reduce this. For example, a planararray in an x−y plane cannot measure wavenumber in the z direction.Because of this, its search grid may be restricted to the hemisphere,0≦θ≦π/2, 0≦φ≦2π, as the lower hemisphere is identical due to theambiguity and provides no additional information.

One approach is to design a set of tapers at each grid location,(θ_(o)±Δθ, φ_(o)±Δφ), based on an eigendecomposition of the matrix

$\begin{matrix}\begin{matrix}{{\underset{\_}{R}}_{\theta_{o}\phi_{o}} = {\frac{1}{4\pi}{\int_{\theta_{o} - {\Delta\theta}}^{\theta_{o} + {\Delta\;\theta}}{\int_{\phi_{o} - {\Delta\;\phi}}^{\phi_{o} + {\Delta\;\phi}}{{\underset{\_}{v}\left( {\theta,\phi} \right)}{{\underset{\_}{v}}^{H}\left( {\theta,\phi} \right)}\sin\;\theta{\mathbb{d}\theta}{\mathbb{d}\phi}}}}}} \\{= {{\underset{\_}{Q}}_{\theta_{o},\phi_{o}}{\underset{\_}{\Lambda}}_{\theta_{o},\phi_{o}}{\underset{\_}{Q}}_{\theta_{o},\phi_{o}}^{H}}}\end{matrix} & (7.28)\end{matrix}$

where v(θ,φ) is the array manifold response vector. The multi-tapers areselected as the eigenvectors corresponding to the D largest eigenvaluesR _(θ) _(o) _(, Φ) _(o) .w _(θ) _(o) _(,Φ) _(o) _(,d) =q _(θ) _(o) _(,Φ) _(o) _(d-1) ,d=1, . . .,D  (7.29)

Alternatively, one may design a single set of tapers, w _(o,d), perhapsat broadside for the array, and steer the main response axis (MRA) ofthis fixed set of weights through angle space.w _(θ) _(o) _(,Φ) _(o) _(,d) =w _(o,d)

v(θ_(o),Φ_(o)),d=1, . . . ,D  (7.30)

Discrete Line Component Processing (Harmonic Analysis)

Defining the Multi-Taper Weight Matrix,W (θ_(o),Φ_(o))=[ w _(θ) _(o) _(,Φ) _(o) _(,1) ,w _(θ) _(o) _(,Φ) _(o)_(,2) , . . . ,w _(θ) _(o) _(,Φ) _(o) _(,D)]  (7.31)

the eigencoefficients are computed for each snapshot, x _(m), asy _(m)(θ_(o),Φ_(o))= W ^(H)(θ_(o),Φ_(o)) x _(m)  (7.32)

The eigencoefficients output for a line component at (θ_(o), φ_(o)) isgiven byq _(θ) _(o) _(,Φ) _(o) =W ^(H)(θ_(o),Φ_(o))1  (7.33)

where 1 is the all one's vector. This vector defines the subspaces usedfor computing a detection statistic.

$\begin{matrix}{{P_{\theta_{o},\phi_{o}} = {{{\underset{\_}{q}}_{\theta_{o},\phi_{o}}\left( {{\underset{\_}{q}}_{\theta_{o},\phi_{o}}^{H}{\underset{\_}{q}}_{\theta_{o},\phi_{o}}} \right)}^{- 1}{\underset{\_}{q}}_{\theta_{o},\phi_{o}}^{H}}},{{\underset{\_}{P}}_{\theta_{o},\phi_{o}}^{\bot} = {\underset{\_}{I} - {\underset{\_}{P}}_{\theta_{o},\phi_{o}}}}} & (7.34)\end{matrix}$

The detection statistic may be computed as before

$\begin{matrix}{{F\left( {\theta_{o},\phi_{o}} \right)} = {\frac{\sum\limits_{m = 1}^{M}{{{\underset{\_}{y}}_{m}^{H}\left( {\theta_{o},\phi_{o}} \right)}{\underset{\_}{P}}_{\theta_{o},\phi_{o}}{{\underset{\_}{y}}_{m}\left( {\theta_{o},\phi_{o}} \right)}}}{\sum\limits_{m = 1}^{M}{{{\underset{\_}{y}}_{m}^{H}\left( {\theta_{o},\phi_{o}} \right)}{\underset{\_}{P}}_{\theta_{o},\phi_{o}}^{\bot}{{\underset{\_}{y}}_{m}\left( {\theta_{o},\phi_{o}} \right)}}}\overset{H_{1}}{\underset{H_{0}}{\gtrless}}{\gamma\;{TH}}}} & (7.35)\end{matrix}$

K is the number of detected line components. Each has its parametersestimated, ({circumflex over (σ)}_(k), {circumflex over (Φ)}_(k),â_(k)(m)), and these are used to subtract the line component from thedata. The projection or subtraction methods described above may be used.The variance for each line component may be estimated as

$\begin{matrix}{{\hat{\sigma}}_{k}^{2} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{{\hat{a}}_{k}(m)}}^{2}}}} & (7.36)\end{matrix}$

and may be used in the formation of the covariance matrix. This processmay be iterated to successively process multiple line components in thedata. The final residual snapshot data, x _(b,m), may include thesmooth, continuous background content.

Background/Continuous Spectrum

Once harmonic analysis is complete the residual snapshot data, x _(b,m),may be used to compute a final smooth, continuous background spectrum.The eigencoefficients are computedy _(b,m)(θ_(o),φ_(o))= W ^(H)(θ_(o),φ_(o)) x _(b,m)=((y _(b,m)^((d))(θ_(o),φ_(o))))_(d)  (7.37)

The eigencoefficients are used to produce the individual eigenspectra.{circumflex over (P)} _(b) ^((d))(θ_(o),φ_(o))=|y _(b,m)^((d))(θ_(o),φ_(o))|²  (7.38)

The individual eigenspectra are then linearly combined according to aset of weights

$\begin{matrix}{{{\hat{P}}_{b,x}\left( {\theta_{o},\phi_{o}} \right)} = {\sum\limits_{d = 1}^{D}{{h_{d}\left( {\theta_{o},\phi_{o}} \right)}{{\hat{P}}_{b}^{(d)}\left( {\theta_{o},\phi_{o}} \right)}}}} & (7.39)\end{matrix}$

The weights, h_(d)(k), may be fixed, which may be optimal for anunderlying white spectrum, or may be determined adaptively.

Estimate Sensor Noise in Virtual Region

As described above, analyzing the covariance matrix of isotropic noiseto determine an appropriate subspace for the virtual region for thearray may provide a dimension, N_(vr), and a projection matrix for thatsubspace, P _(vr). The sensor noise may be estimated from the residualsnapshot data, x _(b,m), as

$\begin{matrix}{{\hat{\sigma}}_{w}^{2} = {{\frac{1}{M} \cdot \frac{1}{N_{vr}}}{\sum\limits_{m = 1}^{M}{{\underset{\_}{x}}_{b,m}^{H}{\underset{\_}{P}}_{vr}{\underset{\_}{x}}_{b,m}}}}} & (7.40)\end{matrix}$

Covariance Matrix Estimate

The final estimate of the covariance matrix may be formed from theestimated line components, the covariance from spatial spectrum of theresidual continuous background, and the sensor noise component as

$\begin{matrix}{{\hat{\underset{\_}{R}}}_{CSS} = {{\hat{\underset{\_}{V}}{\underset{\_}{\hat{R}}}_{a}{\hat{\underset{\_}{V}}}^{H}} + {\left( {4\pi} \right)^{- 1}{\int_{0}^{\pi}{\sin\;\theta{\mathbb{d}\theta}{\int_{\pi}^{\pi}{{{\hat{P}}_{b,x}\left( {\theta,\phi} \right)}{\underset{\_}{v}\left( {\theta,\phi} \right)}{{\underset{\_}{v}}^{H}\left( {\theta,\phi} \right)}{\mathbb{d}\phi}}}}}} + {\sigma_{w}^{2}{\underset{\_}{P}}_{vr}}}} & (7.41)\end{matrix}$

where{circumflex over (R)} _(a)=diag({circumflex over (σ)}₁ ²,{circumflexover (σ)}₂ ², . . . ,{circumflex over (σ)}_(K) ²)  (7.42)and{circumflex over (V)}=[ v ({circumflex over (θ)}₁,{circumflex over(Φ)}₁), v ({circumflex over (θ)}₂,{circumflex over (Φ)}₂), . . . , v({circumflex over (θ)}_(K),{circumflex over (Φ)}_(K))]  (7.43)

Thus, the CSS techniques disclosed herein are applicable to arrays witharbitrary geometry. In a particular embodiment, the sensor noisecomponent in virtual space may be estimated and accounted for. A methodfor doing so based on the covariance matrix of 3D isotropic noise isexplained above. Uniform circular arrays maintain a certain structure,and are not strictly speaking arbitrary, but do result in a covariancematrix which is not Toeplitz, and may be addressed using the samemethods. Normalized SINR loss performance and average MTSE spectra havebeen assessed via simulation. The simulations indicate that performanceis generally very close to optimal with little snapshot support, withsome losses encountered near line components due to mismatch, and forweak INR line components at very low snapshot support due toinconsistent detection within harmonic analysis.

FIG. 6 is a block diagram of a computing environment 600 including ageneral purpose computing device 610 operable to support embodiments ofcomputer-implemented methods and computer-executable programinstructions according to the present disclosure. For example, thecomputing device 610 may include or be included within an arrayprocessor, such as the array processor 102 of FIG. 1. The computingdevice 610 typically includes at least one processor 620. Within thecomputing device 610, the processor 620 communicates with a systemmemory 630, one or more storage devices 640, one or more input/outputinterfaces 650, and one or more communication interfaces 660.

The system memory 630 may include volatile memory devices (e.g., randomaccess memory (RAM) devices), nonvolatile memory devices (e.g.,read-only memory (ROM) devices, programmable read-only memory, and flashmemory), or both. The system memory 630 typically includes an operatingsystem 632, which may include a basic/input output system for bootingthe computing device 610 as well as a full operating system to enablethe computing device 610 to interact with users, other programs, andother devices. The system memory 630 also typically includes one or moreapplication programs 634 and program data 636. For example, theapplication programs 634 and program data 636 may enable the computingdevice 610 to perform one or more aspects of signal processing, asdescribed above. To illustrate, the application programs 634 may includeinstructions that, when executed by the processor 620, cause theprocessor 620 to receive sensed data from sensors of a sensor array(such as the sensor array 110 of FIG. 1. Data from each sensor of thesensor array may be descriptive of waveform phenomena detected at thesensor. The application programs 634 may also include instructions thatcause the processor 620 to determine an estimated spatial spectrum ofthe waveform phenomena based at least partially on the sensed data. Theapplication programs 634 may further include instructions that cause theprocessor 620 to determine an estimated covariance matrix of thewaveform phenomena based on the estimated spatial spectrum. Theapplication programs 634 may also include instructions cause theprocessor 620 to determine adaptive beamforming weights using theestimated covariance matrix.

The processor 620 may also communicate with one or more storage devices640. For example, the one or more storage devices 640 may includenonvolatile storage devices, such as magnetic disks, optical disks, orflash memory devices. The storage devices 640 may include both removableand non-removable memory devices. The storage devices 640 may beconfigured to store an operating system, applications, and program data.In a particular embodiment, the system memory 630, the storage devices640, or both, include tangible, non-transitory computer-readable media.

The processor 620 may also communicate with one or more input/outputinterfaces 650 that enable the computing device 610 to communicate withone or more input/output devices 670 to facilitate user interaction. Theinput/output interfaces 650 may include serial interfaces (e.g.,universal serial bus (USB) interfaces or IEEE 1394 interfaces), parallelinterfaces, display adapters, audio adapters, and other interfaces. Theinput/output devices 670 may include keyboards, pointing devices,displays, speakers, microphones, touch screens, and other devices.

The processor 620 may communicate with other computer systems 680 viathe one or more network interfaces 660. The one or more networkinterfaces 660 may include wired Ethernet interfaces, IEEE 802.01wireless interfaces, Bluetooth communication interfaces, or othernetwork interfaces. The other computer systems 680 may include hostcomputers, servers, workstations, and other computing devices, such asthe other system components 106 of FIG. 1.

The illustrations of the embodiments described herein are intended toprovide a general understanding of the structure of the variousembodiments. The illustrations are not intended to serve as a completedescription of all of the elements and features of apparatus and systemsthat utilize the structures or methods described herein. Many otherembodiments may be apparent to those of skill in the art upon reviewingthe disclosure. Other embodiments may be utilized and derived from thedisclosure, such that structural and logical substitutions and changesmay be made without departing from the scope of the disclosure. Forexample, method steps may be performed in a different order than isshown in the illustrations or one or more method steps may be omitted.Accordingly, the disclosure and the figures are to be regarded asillustrative rather than restrictive.

Moreover, although specific embodiments have been illustrated anddescribed herein, it should be appreciated that any subsequentarrangement designed to achieve the same or similar results may besubstituted for the specific embodiments shown. This disclosure isintended to cover any and all subsequent adaptations or variations ofvarious embodiments. Combinations of the above embodiments and otherembodiments not specifically described herein will be apparent to thoseof skill in the art upon reviewing the description.

In the foregoing Detailed Description, various features may be groupedtogether or described in a single embodiment for the purpose ofstreamlining the disclosure. This disclosure is not to be interpreted asreflecting an intention that the claimed embodiments require morefeatures than are expressly recited in each claim. Rather, as thefollowing claims reflect, the claimed subject matter may be directed toless than all of the features of any of the disclosed embodiments.

What is claimed is:
 1. A computer-implemented method comprising:transforming, by each sensor of a sensor array, electromagnetic energyinto sensed data descriptive of a waveform phenomenon from the sensorsof the sensor array, wherein data from each sensor is descriptive of atleast a portion of the waveform phenomenon, and wherein the sensed dataincludes a discrete component and a continuous component; receiving at aprocessor, the sensed data from the sensors of the sensor array;applying, at the processor, multi-taper spectral estimation and harmonicanalysis to the sensed data to determine an estimated spatial spectrum,wherein applying the harmonic analysis comprises detecting andsubtracting the discrete component from the sensed data; determining, atthe processor, an estimated covariance matrix corresponding to thewaveform phenomenon based on the estimated spatial spectrum;determining, at the processor, adaptive beamforming weights using theestimated covariance matrix; generating, at the processor, a time-domainsignal descriptive of the waveform phenomenon by applying the adaptivebeamforming weights to the sensed data; and sending the time-domainsignal to a computing device.
 2. The method of claim 1, wherein thesensed data corresponds to a single snapshot from the sensor array. 3.The method of claim 1, wherein the harmonic analysis comprises detectingand subtracting a plurality of discrete components from the sensed dataiteratively.
 4. The method of claim 1, wherein a source of the waveformphenomenon is treated as stationary for purposes of determining theestimated covariance matrix.
 5. The method of claim 1, furthercomprising subtracting one or more estimated point source interferencecomponents from the sensed data before determining the estimatedcovariance matrix.
 6. The method of claim 1, wherein the estimatedcovariance matrix is determined based on information regarding locationsof the sensors with respect to one another.
 7. The method of claim 1,wherein the estimated spatial spectrum comprises a frequency-wavenumberspectrum of the waveform phenomenon.
 8. The method of claim 1, whereinthe estimated spatial spectrum comprises a direction-of-arrival spectrumof the waveform phenomenon.
 9. The method of claim 1, wherein the senseddata includes first data acquired by a first sensor of the sensor arrayand second data acquired by a second sensor of the sensor array, andwherein the first data is acquired by the first sensor and the seconddata is acquired by the second sensor simultaneously.
 10. The method ofclaim 1, wherein applying the multi-taper spectral estimation includesusing a Fast Fourier Transform process with zero-padding.
 11. The methodof claim 1, wherein a number of tapers used for the multi-taper spectralestimation is less than a number of the sensors of the sensor array, andwherein zero-padding is used to increase sampling density of theestimated spatial spectrum.
 12. The method of claim 1, wherein thediscrete component is subtracted from the sensed data using a coherentsubtraction method.
 13. The method of claim 1, wherein the harmonicanalysis comprises: detecting the discrete component of the waveformphenomenon using outputs of the multi-taper spectral estimation; anddetermining estimated discrete component parameters corresponding toeach detected discrete component.
 14. The method of claim 13, furthercomprising determining an estimated array error vector corresponding toeach detected discrete component.
 15. The method of claim 14, whereinthe estimated array error vector corresponds to a linear minimum meansquare error estimate of non-ideal response error in the sensed data.16. The method of claim 14, further comprising subtracting acontribution of each detected discrete component from the sensed data togenerate a residual continuous component.
 17. The method of claim 16,further comprising determining an estimated continuous componentcovariance matrix based on the residual continuous component.
 18. Themethod of claim 17, wherein the estimated covariance matrix isdetermined as a numeric combination of the estimated continuouscomponent covariance matrix and an estimated discrete componentcovariance for each detected discrete component.
 19. The method of claim18, wherein the estimated discrete component covariance for a particulardetected discrete component is determined based on the estimated arrayerror vector corresponding to the particular detected discrete componentand based on the estimated discrete component parameters of theparticular detected discrete component.
 20. The method of claim 1,wherein the waveform phenomenon includes a signal of interest, andwherein the sensed data is also descriptive of one or more interferencesignals from one or more interferers.
 21. The method of claim 20,wherein the estimated covariance matrix converges with fewer snapshotsthan interferers.
 22. The method of claim 1, wherein the estimatedcovariance matrix is determined as a combination of a first portion anda second portion, wherein the first portion corresponds to one or morephysically propagating waves of the waveform phenomenon, and wherein thesecond portion is associated with sensor noise.
 23. The method of claim22, further comprising determining an estimate of the sensor noise basedat least partially on the second portion.
 24. The method of claim 1,further comprising: physically moving the sensor array whiletransforming, by the sensors, the electromagnetic energy into the senseddata.
 25. The method of claim 1, further comprising: displaying avisualization based on a spectrum of the time-domain signal.
 26. Asystem comprising: a sensor array including multiple sensors configuredto transform electromagnetic energy into sensed data; and an arrayprocessor coupled to the sensor array, wherein the array processor isconfigured to: receive the sensed data descriptive of a waveformphenomenon from the sensors of the sensor array, wherein data from eachsensor is descriptive of at least a portion of the waveform phenomenon,and wherein the sensed data includes a discrete component and acontinuous component; apply multi-taper spectral estimation and harmonicanalysis to the sensed data to determine an estimated spatial spectrum,wherein applying the multi-taper comprises detecting and subtracting thediscrete component from the sensed data; determine an estimatedcovariance matrix corresponding to the waveform phenomenon based on theestimated spatial spectrum; determine adaptive beamforming weights usingthe estimated covariance matrix; generate a time-domain signaldescriptive of the waveform phenomenon by applying the adaptivebeamforming weights to the sensed data; and send the time-domain signalto a computing device.
 27. The system of claim 26, wherein the harmonicanalysis comprises detecting and subtracting a plurality of discretecomponents from the sensed data iteratively.
 28. The system of claim 26,wherein the sensor array is arranged with non-uniform spacing betweenthe multiple sensors.
 29. A non-transitory computer-readable mediumcomprising instructions that, when executed by a processor, cause theprocessor to: receive sensed data descriptive of a waveform phenomenonfrom sensors of a sensor array, wherein data from each sensor isdescriptive of at least a portion of the waveform phenomenon, andwherein the sensed data includes a discrete component and a continuouscomponent; apply multi-taper spectral estimation and harmonic analysisto the sensed data to determine an estimated spatial spectrum, whereinapplying the harmonic analysis comprises detecting and subtracting thediscrete component from the sensed data; determine an estimatedcovariance matrix corresponding to the waveform phenomenon based on theestimated spatial spectrum; determine adaptive beamforming weights usingthe estimated covariance matrix; generate a time-domain signaldescriptive of the waveform phenomenon by applying the adaptivebeamforming weights to the sensed data; and send the time-domain signalto a computing device.
 30. The non-transitory computer-readable mediumof claim 29, wherein the sensed data corresponds to a single snapshotfrom the sensor array.